### Replication files for
### "Local Immigration Policies Shape Immigrants' Relocation Preferences Regardless of Immigration Status"
### Hans Lueders
### Political Behavior
### August 2025




# Set-up ------------------------------------------------------------------

### packages
library(ggplot2)
library(foreign)
require(lmtest)
library(stargazer)
library(xtable)
library(maptools)
library(gridExtra)
library(ggmap)
library(rgdal)
library(sp)
library(gridExtra)
library(grid)
library(lattice)
library(Matching)
library(ggpubr)
library(geosphere)
library(spatstat)
library(raster)
library(lfe)
library(AER)
library(dplyr)
library(corrplot)
library(ggrepel)
library(rbounds)
library(mgcv)
library(mediation)
library(sensemakr)
library(cjoint)
library(emmeans)
library(cregg)
library(mapproj)
library(readxl)
library(sf)
library(Hmisc)
library(estimatr)



### empty environment
rm(list=ls())


### set working directory
path <- "" # define path here
setwd(path)






# Load and prepare datasets -----------------------------------------------

### Load data
load(paste(path, "LocalImmigrationPolicies_replicationdata.Rdata", sep=""))



# Figure 1: Immigration-related concerns by immigration status ------------

### Create output container
out <- data.frame(matrix(nrow = 12, ncol=4))
colnames(out) <- c("variable", "group", "mean", "se")
out$variable <- factor(rep(1:4, each=3), levels = 1:4,
                       labels = c("Worried\nabout deportation", 
                                  "Even legal immigrants\nshould worry about deportation",
                                  "Unauthorized\nfamily member",
                                  "Know someone\nwho was deported"))
out$group <- factor(rep(1:3,4), levels = 1:3,
                    labels = c("Naturalized", "LPR", "Likely unauthorized"))


### get values
out[1:3,3:4] <- summary(lm_robust(cit_worries_deportation_dich ~ cit_citizenUS_alt - 1, data=d))$coef[2:4,1:2]
out[4:6,3:4] <- summary(lm_robust(bl_agree_deportation_dich ~ cit_citizenUS_alt - 1, data=d))$coef[2:4,1:2]
out[7:9,3:4] <- summary(lm_robust(cit_family_unauthorized ~ cit_citizenUS_alt - 1, data=d))$coef[2:4,1:2]
out[10:12,3:4] <- summary(lm_robust(cit_know_deported ~ cit_citizenUS_alt - 1, data=d))$coef[2:4,1:2]


### Figure
p <- ggplot() +
  geom_bar(data=out, 
           aes(x=group, y=mean*100, fill=group), 
           stat="identity", alpha=.8) +
  geom_errorbar(data=out, 
                aes(x=group, ymin=(mean-1.96*se)*100, ymax=(mean+1.96*se)*100, col=group),
                width=0.15, size=1) +
  scale_y_continuous(limits = c(0,86), breaks=seq(0,80,20)) +
  scale_color_manual(values = c("burlywood3", "peachpuff4", "gray20"), name="") +
  scale_fill_manual(values = c("burlywood2", "peachpuff3", "gray30"), name="") + 
  facet_wrap(~variable) +
  ylab("Share of respondents") +
  theme_bw() + 
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=10),
        legend.position = "none")  
p
ggsave("Figure1_Immiworries_unweighted.png", p, width=10, height=5)




# Figure 2: Results of the push experiment --------------------------------


### create output container
out <- data.frame(matrix(nrow=8*5, ncol=5))
colnames(out) <- c("status","group", "condition", "b", "se")
out$status <- rep(1:5, each=8)
out$status <- factor(out$status, levels=1:5,
                     labels = c("all", "all (complete)", "naturalized", "LPR", "likely unauthorized"))
out$group <- rep(c(1,1,1,1,2,2,2,2),5)
out$group <- factor(out$group, levels=1:2,
                    labels = c("positive\nchange", "negative\nchange"))
out$condition <- rep(c(1,2,3,4,1,2,3,4),5)
out$condition <- factor(out$condition, levels=1:4,
                        labels = c("Economic opportunities", "Social network", "Deportation risk", "Healthcare access"))


### get values
out[1:4,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d.full))$coef[1:4,1:2]
out[5:8,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d.full))$coef[1:4,1:2]

out[9:12,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d))$coef[1:4,1:2]
out[13:16,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d))$coef[1:4,1:2]  

out[17:20,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, 
                                    subset=d$cit_citizenUS_alt == "naturalized US citizen"))$coef[1:4,1:2]
out[21:24,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, 
                                    subset=d$cit_citizenUS_alt == "naturalized US citizen"))$coef[1:4,1:2]

out[25:28,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, 
                                    subset=d$cit_citizenUS_alt == "lawful permanent resident"))$coef[1:4,1:2]
out[29:32,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, 
                                    subset=d$cit_citizenUS_alt == "lawful permanent resident"))$coef[1:4,1:2]  

out[33:36,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, 
                                    subset=d$cit_citizenUS_alt == "likely unauthorized immigrant"))$coef[1:4,1:2]
out[37:40,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, 
                                    subset=d$cit_citizenUS_alt == "likely unauthorized immigrant"))$coef[1:4,1:2] 


### plot
p <- ggplot() + 
  geom_bar(data=out, aes(x=group, y=b, fill=status), 
           stat="identity", position = position_dodge(), alpha=.9) +
  geom_errorbar(data=out, aes(x=group, ymin=b-1.96*se, ymax=b+1.96*se,col=status),
                position = position_dodge(width=0.9), width=0.25, size=1) +
  scale_fill_manual(values = c("gray40", "gray60", "burlywood1","burlywood3","burlywood4"), name="") +
  scale_color_manual(values = rep("gray30",5), name="") +
  geom_hline(yintercept = 0, lty=2) +
  xlab("") + 
  #ylab("Change in Pr(move in next five years)\n(neg=less likely; pos=more likely)") +
  ylab("Change in Pr(move)") +
  facet_wrap(~condition) +  
  theme_bw() +
  theme(strip.text = element_text(size=14,  face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12)) 
p
ggsave("Figure2_Push_results_main.png", p, width=8, height=6)




# Figure 3: Immigration-related concerns moderate the impact of lo --------

### create output container
out <- data.frame(matrix(nrow=16, ncol=5))
colnames(out) <- c("direction", "moderator", "level", "b", "se")
out$direction <- factor(rep(rep(1:2,each=2),4), levels = 1:2,
                        labels = c("positive\nchange", "negative\nchange"))
out$moderator <- factor(rep(1:4,each=4), levels = 1:4,
                        labels = c("Worried\nabout deportation", 
                                   "Even legal immigrants\nshould worry about deportation",
                                   "Unauthorized\nfamily member","Know someone\nwho was deported"))
out$level <- factor(rep(1:2,8), levels = 1:2,
                    labels = c("no", "yes"))

### get values
out[1,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$cit_worries_deportation_dich == 0))$coef[3,1:2]
out[2,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$cit_worries_deportation_dich == 1))$coef[3,1:2]

out[3,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$cit_worries_deportation_dich == 0))$coef[3,1:2]
out[4,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$cit_worries_deportation_dich == 1))$coef[3,1:2]

out[5,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$bl_agree_deportation_dich == 0))$coef[3,1:2]
out[6,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$bl_agree_deportation_dich == 1))$coef[3,1:2]

out[7,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$bl_agree_deportation_dich == 0))$coef[3,1:2]
out[8,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$bl_agree_deportation_dich == 1))$coef[3,1:2]

out[9,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$cit_family_unauthorized == 0))$coef[3,1:2]
out[10,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$cit_family_unauthorized == 1))$coef[3,1:2]

out[11,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$cit_family_unauthorized == 0))$coef[3,1:2]
out[12,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$cit_family_unauthorized == 1))$coef[3,1:2]

out[13,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$cit_know_deported == 0))$coef[3,1:2]
out[14,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$cit_know_deported == 1))$coef[3,1:2]

out[15,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$cit_know_deported == 0))$coef[3,1:2]
out[16,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$cit_know_deported == 1))$coef[3,1:2]


### plot
p <- ggplot() + 
  geom_bar(data=out, 
           aes(x=as.numeric(direction), y=b, fill=level), 
           stat="identity", position = position_dodge(), alpha=.9) +
  geom_errorbar(data=out, 
                aes(x=as.numeric(direction), ymin=b-1.96*se, ymax=b+1.96*se,col=level),
                position = position_dodge(width=0.9), width=0.25, size=1) +
  scale_x_continuous(breaks=c(1,2), labels = c("positive\nchange", "negative\nchange")) +
  scale_fill_manual(values = c("burlywood2","burlywood4"), name="") +
  scale_color_manual(values = c("gray30","gray30"), name="") +
  geom_hline(yintercept = 0, lty=2) +
  xlab("") + 
  ylab("Change in Pr(move)") +
  #ylab("Change in Pr(move in next five years)\n(neg=less likely; pos=more likely)") +
  facet_wrap(~moderator) +  
  theme_bw() +
  theme(strip.text = element_text(size=14,  face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12)) 
p
ggsave("Figure3_Push_results_HTE.png", p, width=8, height=6)




# Figure 4: Interaction with state immigration environment ----------------

### create output container
out <- data.frame(matrix(nrow=8, ncol=5))
colnames(out) <- c("direction", "moderator", "level", "b", "se")
out$direction <- factor(rep(rep(1:2,each=2),2), levels = 1:2,
                        labels = c("positive\nchange", "negative\nchange"))
out$moderator <- factor(rep(1:2,each=4), levels = 1:2,
                        labels = c("Democratic\nGovernor", "Democratic\nState Legislature"))
out$level <- factor(rep(1:2,4), levels = 1:2,
                    labels = c("no", "yes"))


### get values
out[1,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$DemGovernor == 0))$coef[3,1:2]
out[2,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$DemGovernor == 1))$coef[3,1:2]

out[3,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$DemGovernor == 0))$coef[3,1:2]
out[4,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$DemGovernor == 1))$coef[3,1:2]

out[5,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$DemLegislature == 0))$coef[3,1:2]
out[6,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=d, subset=d$DemLegislature == 1))$coef[3,1:2]

out[7,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$DemLegislature == 0))$coef[3,1:2]
out[8,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=d, subset=d$DemLegislature == 1))$coef[3,1:2]


### Plot
p <- ggplot() + 
  geom_bar(data=out, 
           aes(x=as.numeric(direction), y=b, fill=level), 
           stat="identity", position = position_dodge(), alpha=.9) +
  geom_errorbar(data=out, 
                aes(x=as.numeric(direction), ymin=b-1.96*se, ymax=b+1.96*se,col=level),
                position = position_dodge(width=0.9), width=0.25, size=1) +
  scale_x_continuous(breaks=c(1,2), labels = c("positive\nchange", "negative\nchange")) +
  scale_fill_manual(values = c("burlywood2","burlywood4"), name="") +
  scale_color_manual(values = c("gray30","gray30"), name="") +
  geom_hline(yintercept = 0, lty=2) +
  xlab("") + 
  ylab("Change in Pr(move)") +
  #ylab("Change in Pr(move in next five years)\n(neg=less likely; pos=more likely)") +
  facet_wrap(~moderator) +  
  theme_bw() +
  theme(strip.text = element_text(size=14,  face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12)) 
p
ggsave("Figure4_Push_results_state.png", p, width=7, height=4)





# Figure 5: Determinants of destination preferences -----------------------

### Plot 1: all respondents
# prepare output container
out <- data.frame(matrix(nrow=35*2, ncol=6))
colnames(out) <- c("sample","dimension", "attribute", "b", "lb", "ub")
out$sample <- factor(rep(1:2,each=35), levels = 1:2,
                     labels = c("all", "all (complete)"))
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- factor(rep(35:1,2), levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients: all respondents
reg <- mm(data = d.full.conjoint, formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# get coefficients
out[2:5,4:6] <- reg[1:4,c(5,9,10)]
out[7:10,4:6] <- reg[5:8,c(5,9,10)]
out[12:15,4:6] <- reg[9:12,c(5,9,10)]
out[17:20,4:6] <- reg[13:16,c(5,9,10)]
out[22:25,4:6] <- reg[17:20,c(5,9,10)]
out[27:30,4:6] <- reg[21:24,c(5,9,10)]
out[32:35,4:6] <- reg[25:28,c(5,9,10)]

# get coefficients: all respondents, completes
reg <- mm(data = d.conjoint, formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# get coefficients
out[37:40,4:6] <- reg[1:4,c(5,9,10)]
out[42:45,4:6] <- reg[5:8,c(5,9,10)]
out[47:50,4:6] <- reg[9:12,c(5,9,10)]
out[52:55,4:6] <- reg[13:16,c(5,9,10)]
out[57:60,4:6] <- reg[17:20,c(5,9,10)]
out[62:65,4:6] <- reg[21:24,c(5,9,10)]
out[67:70,4:6] <- reg[25:28,c(5,9,10)]


# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             position=position_dodge(width=0.8)) +
  #geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
  #           shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +  
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_shape_manual(values = c(17,19), name="") +  
  scale_color_manual(values = c("chocolate3","chocolate4"), name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("Marginal Mean") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


### Plot 2: by citizenship
# Regression
mm_by <- cj(d.conjoint, 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +
  scale_linetype_manual(values = c(2,1), name="") +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



### Plot 3: Difference in marginal means between naturalized citizens vs. likely unauthorized immigrants
# dichotomize moderator
d.conjoint$citDich <- NA
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "naturalized US citizen"] <- 0
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "likely unauthorized immigrant"] <- 1
d.conjoint$citDich <- factor(d.conjoint$citDich, levels=0:1,
                             labels = c("naturalized", "likely unauthorized"))

# Regression
mm_diff <- cj(d.conjoint, 
              E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_differences",
              by=~citDich)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- mm_diff[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_diff[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_diff[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_diff[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_diff[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_diff[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_diff[25:28,c(6,10,11)]

# define sample
out$sample <- "difference"

# define comparison
out$comparison <- "Likely unauth. - naturalized"

# Plot
p3 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
             #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("chocolate4"), name="") +
  scale_y_continuous(limits = c(-0.15,0.15)) +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Difference in Marginal Means") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



### combining all three plots
p <- ggarrange(p1,p2,p3,widths=c(0.45,0.25,0.25), nrow = 1)
p
ggsave("Figure5_Pull_results_main.png", p, width=13, height=7)





# Figure 6: Determinants of destination preferences, by immigratio --------

### Have unauthorized family member
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_family_unauthorized)

# prepare output container
out1 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out1) <- c("moderator","sample","dimension", "attribute", "b", "lb", "ub")
out1$moderator <- 3 # "Unauthorized\nfamily member"
out1$sample <- rep(1:2,each=35)
out1$sample <- factor(out1$sample, levels = 1:2,
                      labels = c("no", "yes"))
out1$dimension <- rep(rep(1:7, each=5),2)
out1$dimension <- factor(out1$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out1$attribute <- rep(35:1,2)
out1$attribute <- factor(out1$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out1[2:5,5:7] <- reg[1:4,c(6,10,11)]
out1[7:10,5:7] <- reg[5:8,c(6,10,11)]
out1[12:15,5:7] <- reg[9:12,c(6,10,11)]
out1[17:20,5:7] <- reg[13:16,c(6,10,11)]
out1[22:25,5:7] <- reg[17:20,c(6,10,11)]
out1[27:30,5:7] <- reg[21:24,c(6,10,11)]
out1[32:35,5:7] <- reg[25:28,c(6,10,11)]

out1[37:40,5:7] <- reg[29:32,c(6,10,11)]
out1[42:45,5:7] <- reg[33:36,c(6,10,11)]
out1[47:50,5:7] <- reg[37:40,c(6,10,11)]
out1[52:55,5:7] <- reg[41:44,c(6,10,11)]
out1[57:60,5:7] <- reg[45:48,c(6,10,11)]
out1[62:65,5:7] <- reg[49:52,c(6,10,11)]
out1[67:70,5:7] <- reg[53:56,c(6,10,11)]



### Know someone who was deported
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_know_deported)

# prepare output container
out2 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out2) <- c("moderator","sample","dimension", "attribute", "b", "lb", "ub")
out2$moderator <- 4 #"Knows someone\nwho was deported"
out2$sample <- rep(1:2,each=35)
out2$sample <- factor(out2$sample, levels = 1:2,
                      labels = c("no", "yes"))
out2$dimension <- rep(rep(1:7, each=5),2)
out2$dimension <- factor(out2$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out2$attribute <- rep(35:1,2)
out2$attribute <- factor(out2$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out2[2:5,5:7] <- reg[1:4,c(6,10,11)]
out2[7:10,5:7] <- reg[5:8,c(6,10,11)]
out2[12:15,5:7] <- reg[9:12,c(6,10,11)]
out2[17:20,5:7] <- reg[13:16,c(6,10,11)]
out2[22:25,5:7] <- reg[17:20,c(6,10,11)]
out2[27:30,5:7] <- reg[21:24,c(6,10,11)]
out2[32:35,5:7] <- reg[25:28,c(6,10,11)]

out2[37:40,5:7] <- reg[29:32,c(6,10,11)]
out2[42:45,5:7] <- reg[33:36,c(6,10,11)]
out2[47:50,5:7] <- reg[37:40,c(6,10,11)]
out2[52:55,5:7] <- reg[41:44,c(6,10,11)]
out2[57:60,5:7] <- reg[45:48,c(6,10,11)]
out2[62:65,5:7] <- reg[49:52,c(6,10,11)]
out2[67:70,5:7] <- reg[53:56,c(6,10,11)]



### Even legal immigrants should be worried about deportation
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~bl_agree_deportation_dich)

# prepare output container
out3 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out3) <- c("moderator","sample","dimension", "attribute", "b", "lb", "ub")
out3$moderator <- 2 #"Even legal immigrants should worry about deportation"
out3$sample <- rep(1:2,each=35)
out3$sample <- factor(out3$sample, levels = 1:2,
                      labels = c("no", "yes"))
out3$dimension <- rep(rep(1:7, each=5),2)
out3$dimension <- factor(out3$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out3$attribute <- rep(35:1,2)
out3$attribute <- factor(out3$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out3[2:5,5:7] <- reg[1:4,c(6,10,11)]
out3[7:10,5:7] <- reg[5:8,c(6,10,11)]
out3[12:15,5:7] <- reg[9:12,c(6,10,11)]
out3[17:20,5:7] <- reg[13:16,c(6,10,11)]
out3[22:25,5:7] <- reg[17:20,c(6,10,11)]
out3[27:30,5:7] <- reg[21:24,c(6,10,11)]
out3[32:35,5:7] <- reg[25:28,c(6,10,11)]

out3[37:40,5:7] <- reg[29:32,c(6,10,11)]
out3[42:45,5:7] <- reg[33:36,c(6,10,11)]
out3[47:50,5:7] <- reg[37:40,c(6,10,11)]
out3[52:55,5:7] <- reg[41:44,c(6,10,11)]
out3[57:60,5:7] <- reg[45:48,c(6,10,11)]
out3[62:65,5:7] <- reg[49:52,c(6,10,11)]
out3[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Worried about deportation
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_worries_deportation_dich)

# prepare output container
out4 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out4) <- c("moderator","sample","dimension", "attribute", "b", "lb", "ub")
out4$moderator <- 1 #"Worried\nabout deportation"
out4$sample <- rep(1:2,each=35)
out4$sample <- factor(out4$sample, levels = 1:2,
                      labels = c("no", "yes"))
out4$dimension <- rep(rep(1:7, each=5),2)
out4$dimension <- factor(out4$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out4$attribute <- rep(35:1,2)
out4$attribute <- factor(out4$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out4[2:5,5:7] <- reg[1:4,c(6,10,11)]
out4[7:10,5:7] <- reg[5:8,c(6,10,11)]
out4[12:15,5:7] <- reg[9:12,c(6,10,11)]
out4[17:20,5:7] <- reg[13:16,c(6,10,11)]
out4[22:25,5:7] <- reg[17:20,c(6,10,11)]
out4[27:30,5:7] <- reg[21:24,c(6,10,11)]
out4[32:35,5:7] <- reg[25:28,c(6,10,11)]

out4[37:40,5:7] <- reg[29:32,c(6,10,11)]
out4[42:45,5:7] <- reg[33:36,c(6,10,11)]
out4[47:50,5:7] <- reg[37:40,c(6,10,11)]
out4[52:55,5:7] <- reg[41:44,c(6,10,11)]
out4[57:60,5:7] <- reg[45:48,c(6,10,11)]
out4[62:65,5:7] <- reg[49:52,c(6,10,11)]
out4[67:70,5:7] <- reg[53:56,c(6,10,11)]



### combine all three
out <- rbind(out1,out2)
out <- rbind(out,out3)
out <- rbind(out,out4)


### factorize moderators
out$moderator <- factor(out$moderator, levels = 1:4,
                        labels = c("Worried\nabout deportation",
                                   "Even legal immigrants\nshould worry\nabout deportation",
                                   "Unauthorized\nfamily member",
                                   "Knows someone\nwho was deported"))




### plot
p <- ggplot() + 
  geom_hline(yintercept = 50, lty=2) +
  geom_point(data=out[out$dimension %in% c("Deportation risk"),], 
             aes(x=attribute, y=b*100, shape=sample, col=sample), size=2.5,
             position = position_dodge(width=0.5)) +
  geom_errorbar(data=out[out$dimension %in% c("Deportation risk"),], 
                aes(x=attribute, ymin=lb*100, ymax=ub*100, col=sample), width=0,
                position = position_dodge(width=0.5),size=.75) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  scale_y_continuous(breaks = seq(20,80,10)) +
  facet_wrap(~moderator, nrow=1, scales="free_x") +
  coord_flip() + 
  xlab("") +
  ylab("Marginal Mean") + 
  theme_bw() +  
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_text(size=14, face="bold"),
        axis.text.y = element_text(face=rep(c(rep("plain",4), "bold"),2), size=12),
        axis.text.x = element_text(size=12),
        axis.ticks.y = element_blank(),
        strip.text = element_text(size=14, face="bold"))
p
ggsave("Figure6_Pull_results_HTE.png", p, width=13, height=4)



# Figure 7: Effects of immigration policies on other outcomes -------------

### Fear discrimination
# regression
reg <- cj(d.conjoint, 
          E2_discrimination ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out1 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out1) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out1$outcome <- "Fear\ndiscrimination"
out1$sample <- rep(1:2,each=35)
out1$sample <- factor(out1$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out1$dimension <- rep(rep(1:7, each=5),2)
out1$dimension <- factor(out1$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out1$attribute <- rep(35:1,2)
out1$attribute <- factor(out1$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out1[2:5,5:7] <- reg[1:4,c(6,10,11)]
out1[7:10,5:7] <- reg[5:8,c(6,10,11)]
out1[12:15,5:7] <- reg[9:12,c(6,10,11)]
out1[17:20,5:7] <- reg[13:16,c(6,10,11)]
out1[22:25,5:7] <- reg[17:20,c(6,10,11)]
out1[27:30,5:7] <- reg[21:24,c(6,10,11)]
out1[32:35,5:7] <- reg[25:28,c(6,10,11)]

out1[37:40,5:7] <- reg[29:32,c(6,10,11)]
out1[42:45,5:7] <- reg[33:36,c(6,10,11)]
out1[47:50,5:7] <- reg[37:40,c(6,10,11)]
out1[52:55,5:7] <- reg[41:44,c(6,10,11)]
out1[57:60,5:7] <- reg[45:48,c(6,10,11)]
out1[62:65,5:7] <- reg[49:52,c(6,10,11)]
out1[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Be safe
# regression
reg <- cj(d.conjoint, 
          E2_besafe ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out2 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out2) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out2$outcome <- "Be safe"
out2$sample <- rep(1:2,each=35)
out2$sample <- factor(out2$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out2$dimension <- rep(rep(1:7, each=5),2)
out2$dimension <- factor(out2$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out2$attribute <- rep(35:1,2)
out2$attribute <- factor(out2$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out2[2:5,5:7] <- reg[1:4,c(6,10,11)]
out2[7:10,5:7] <- reg[5:8,c(6,10,11)]
out2[12:15,5:7] <- reg[9:12,c(6,10,11)]
out2[17:20,5:7] <- reg[13:16,c(6,10,11)]
out2[22:25,5:7] <- reg[17:20,c(6,10,11)]
out2[27:30,5:7] <- reg[21:24,c(6,10,11)]
out2[32:35,5:7] <- reg[25:28,c(6,10,11)]

out2[37:40,5:7] <- reg[29:32,c(6,10,11)]
out2[42:45,5:7] <- reg[33:36,c(6,10,11)]
out2[47:50,5:7] <- reg[37:40,c(6,10,11)]
out2[52:55,5:7] <- reg[41:44,c(6,10,11)]
out2[57:60,5:7] <- reg[45:48,c(6,10,11)]
out2[62:65,5:7] <- reg[49:52,c(6,10,11)]
out2[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Local government protects immigrants
# regression
reg <- cj(d.conjoint, 
          E2_govprotect ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out3 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out3) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out3$outcome <- "Gov't protects\nimmigrants"
out3$sample <- rep(1:2,each=35)
out3$sample <- factor(out3$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out3$dimension <- rep(rep(1:7, each=5),2)
out3$dimension <- factor(out3$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out3$attribute <- rep(35:1,2)
out3$attribute <- factor(out3$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out3[2:5,5:7] <- reg[1:4,c(6,10,11)]
out3[7:10,5:7] <- reg[5:8,c(6,10,11)]
out3[12:15,5:7] <- reg[9:12,c(6,10,11)]
out3[17:20,5:7] <- reg[13:16,c(6,10,11)]
out3[22:25,5:7] <- reg[17:20,c(6,10,11)]
out3[27:30,5:7] <- reg[21:24,c(6,10,11)]
out3[32:35,5:7] <- reg[25:28,c(6,10,11)]

out3[37:40,5:7] <- reg[29:32,c(6,10,11)]
out3[42:45,5:7] <- reg[33:36,c(6,10,11)]
out3[47:50,5:7] <- reg[37:40,c(6,10,11)]
out3[52:55,5:7] <- reg[41:44,c(6,10,11)]
out3[57:60,5:7] <- reg[45:48,c(6,10,11)]
out3[62:65,5:7] <- reg[49:52,c(6,10,11)]
out3[67:70,5:7] <- reg[53:56,c(6,10,11)]

### Government is liberal
# regression
reg <- cj(d.conjoint, 
          E2_govliberal ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out4 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out4) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out4$outcome <- "Gov't is\nliberal"
out4$sample <- rep(1:2,each=35)
out4$sample <- factor(out4$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out4$dimension <- rep(rep(1:7, each=5),2)
out4$dimension <- factor(out4$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out4$attribute <- rep(35:1,2)
out4$attribute <- factor(out4$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out4[2:5,5:7] <- reg[1:4,c(6,10,11)]
out4[7:10,5:7] <- reg[5:8,c(6,10,11)]
out4[12:15,5:7] <- reg[9:12,c(6,10,11)]
out4[17:20,5:7] <- reg[13:16,c(6,10,11)]
out4[22:25,5:7] <- reg[17:20,c(6,10,11)]
out4[27:30,5:7] <- reg[21:24,c(6,10,11)]
out4[32:35,5:7] <- reg[25:28,c(6,10,11)]

out4[37:40,5:7] <- reg[29:32,c(6,10,11)]
out4[42:45,5:7] <- reg[33:36,c(6,10,11)]
out4[47:50,5:7] <- reg[37:40,c(6,10,11)]
out4[52:55,5:7] <- reg[41:44,c(6,10,11)]
out4[57:60,5:7] <- reg[45:48,c(6,10,11)]
out4[62:65,5:7] <- reg[49:52,c(6,10,11)]
out4[67:70,5:7] <- reg[53:56,c(6,10,11)]


### combine all four
out <- rbind(out1,out2)
out <- rbind(out,out3)
out <- rbind(out,out4)


### Plot
p <- ggplot() + 
  geom_hline(yintercept = 50, lty=2) +
  geom_point(data=out[out$dimension %in% c("Deportation risk"),], 
             aes(x=attribute, y=b*100, shape=sample, col=sample), size=2.5,
             position = position_dodge(width=0.5)) +
  geom_errorbar(data=out[out$dimension %in% c("Deportation risk"),], 
                aes(x=attribute, ymin=lb*100, ymax=ub*100, col=sample), width=0,
                position = position_dodge(width=0.5),size=.75) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  scale_y_continuous(breaks = seq(20,80,10)) +
  facet_wrap(~outcome, nrow=1, scales="free_x") +
  coord_flip() + 
  xlab("") +
  ylab("Marginal Mean") + 
  theme_bw() +  
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_text(size=14, face="bold"),
        axis.text.y = element_text(face=rep(c(rep("plain",4), "bold"),2), size=12),
        axis.text.x = element_text(size=12),
        axis.ticks.y = element_blank(),
        strip.text = element_text(size=14, face="bold"))
p
ggsave("Figure7_Pull_results_other.png", p, width=11, height=3.5)






# Figure A2: Map of all respondents ---------------------------------------



### Prepare shapefile (including moving Alaska and Hawaii)
# read U.S. counties moderately-simplified GeoJSON file
us <- read_sf(paste(path, "/US shapefile/cb_2018_us_state_5m.shp", sep=""))

# drop non-states
us <- subset(us, !NAME %in% c("United States Virgin Islands", "Guam", "Commonwealth of the Northern Mariana Islands", "American Samoa", "Puerto Rico"))

# convert it to Albers equal area
us_aea <- st_transform(us, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))

# extract, then rotate, shrink & move alaska (and reset projection)
# need to use state IDs via # https://www.census.gov/geo/reference/ansi_statetables.html
alaska <- us_aea %>% filter(NAME %in% 'Alaska')
alaska_g <- st_geometry(alaska)
alaska_centroid <- st_centroid(st_union(alaska_g))
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
alaska_trans <- (alaska_g - alaska_centroid) * rot(-39 * pi/180) / 2.3 + alaska_centroid + c(1000000, -5000000)
alaska <- alaska %>% st_set_geometry(alaska_trans) %>% st_set_crs(st_crs(us_aea))

# same for Hawaii
hawaii <- us_aea %>% filter(NAME %in% 'Hawaii')
hawaii_g <- st_geometry(hawaii)
hawaii_centroid <- st_centroid(st_union(hawaii_g))
hawaii_trans <- (hawaii_g - hawaii_centroid) * rot(-35 * pi/180) + hawaii_centroid + c(5200000, -1400000)
hawaii <- hawaii %>% st_set_geometry(hawaii_trans) %>% st_set_crs(st_crs(us_aea))

# combining everything
map_final <- us_aea %>%
  filter(!NAME %in% c('Alaska', 'Hawaii')) %>%
  rbind(alaska) %>%
  rbind(hawaii)

### prepare coordinates
sp <- SpatialPoints(d[!is.na(d$demo_zipcode_longitude) & !is.na(d$demo_zipcode_state),
                      c("demo_zipcode_longitude", "demo_zipcode_latitude")], 
                    proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
sp2 <- spTransform(sp, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))

sp2.df <- as.data.frame(sp2)
sp2.df <- subset(sp2.df, demo_zipcode_longitude < 9000000)


### Map
p <- ggplot() + 
  geom_sf(data=map_final) +
  geom_point(data=sp2.df, aes(x=demo_zipcode_longitude, y=demo_zipcode_latitude)) +
  theme_bw() + 
  theme(axis.title=element_blank(), 
        axis.text=element_blank(), 
        axis.ticks=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        legend.position="bottom",
        legend.title = element_text(size=14, face="bold"),
        legend.text = element_text(size=12),
        plot.title = element_text(size=18, face="bold", hjust=.5)) + 
  guides(fill = guide_legend(title.position = "top", title.hjust =0.5))
p
ggsave("FigureA2_RespondentsMap.png",p,width=10, height=8)




# Table A3: Response quality ----------------------------------------------

### Create output container
out <- data.frame(matrix(nrow = 6, ncol=2))
colnames(out) <- c("Criterion", "Share")
out$Criterion <- c("Very fast response", 
                   "Very slow response",
                   "Failed attention check",
                   "Inconsistent answer 1: Age",
                   "Inconsistent answer 2: HH size",
                   "Inconsistent answer 3: Conjoint")

### get values
out$Share[1] <- prop.table(table(d$duration < 0.5*median(d$duration, na.rm=T)))[2]
out$Share[2] <- prop.table(table(d$duration > 2*median(d$duration, na.rm=T)))[2]
out$Share[3] <- prop.table(table(d$attentioncheck))[1]
out$Share[4] <- 0 #prop.table(table(d$cit_age_moveUS > d$demo_age))
out$Share[5] <- prop.table(table(d$demo_nHH_adults == 0))[2]
out$Share[6] <- prop.table(table(d.conjoint$E2_choice_consistent))[1]

### convert to percent
out$Share <- out$Share * 100

### Table
print(xtable(out, caption = "Response quality", align = c("lrc")),
      caption.placement = "top", include.rownames=FALSE)




# Table A4: Sample demographics compared to population benchmarks ---------

### Get values
# Gender
out1 <- cbind(prop.table(table(d$demo_gender[d$cit_dich == "non-citizen"])), prop.table(table(d$demo_gender[d$cit_dich == "naturalized citizen"])),
              c(0.5394,0.4606), c(0.4944,0.5056))
out1 <- data.frame(cbind("Gender",rownames(out1),out1))
colnames(out1) <- c("variable","value","sample_noncitizen", "sample_naturalized", "benchmark_noncitizen", "benchmark_naturalized")

# Age group
out2 <- cbind(prop.table(table(d$demo_agegroupALT[d$cit_dich == "non-citizen"])), prop.table(table(d$demo_agegroupALT[d$cit_dich == "naturalized citizen"])),
              c(0.1575,0.2523,0.2914,0.1773,0.0781,0.0435), c(0.0751,0.1511,0.2124,0.2567,0.1769,0.1279))
out2 <- data.frame(cbind("Age group",rownames(out2),out2))
colnames(out2) <- c("variable","value","sample_noncitizen", "sample_naturalized", "benchmark_noncitizen", "benchmark_naturalized")

# Education
out3 <- cbind(prop.table(table(d$demo_educationALT[d$cit_dich == "non-citizen"])), prop.table(table(d$demo_educationALT[d$cit_dich == "naturalized citizen"])),
              c(0.1449,0.6775,0.0939,0.0247,0.0445,0.0144), c(0.1245,0.5499,0.1498,0.0534,0.0859,0.0365))
out3 <- data.frame(cbind("Education",rownames(out3),out3))
colnames(out3) <- c("variable","value","sample_noncitizen", "sample_naturalized", "benchmark_noncitizen", "benchmark_naturalized")

# Income
out4 <- cbind(prop.table(table(d$demo_income_quintile[d$cit_dich == "non-citizen"])), prop.table(table(d$demo_income_quintile[d$cit_dich == "naturalized citizen"])),
              c(0.0633,0.2468,0.2485,0.2216,0.2198), c(0.0596,0.1756,0.2021,0.2300,0.3326))
out4 <- data.frame(cbind("Household income",rownames(out4),out4))
colnames(out4) <- c("variable","value","sample_noncitizen", "sample_naturalized", "benchmark_noncitizen", "benchmark_naturalized")


# combine all
out <- rbind(out1,out2,out3,out4)
rownames(out) <- NULL
out$sample_noncitizen <- as.numeric(out$sample_noncitizen) * 100
out$sample_naturalized <- as.numeric(out$sample_naturalized) * 100
out$benchmark_noncitizen <- as.numeric(out$benchmark_noncitizen) * 100
out$benchmark_naturalized <- as.numeric(out$benchmark_naturalized) * 100


### Output table
print(xtable(out, caption = "Sample demographics compared to population benchmarks, by citizenship status", align = c("lrrcccc")),
      caption.placement = "top", include.rownames=FALSE)



# Table A5: Comparison of select statistics by immigration status ---------

### function to extract statistics
cit.fun <- function(var,varname){
  # turn variable into numeric
  d$dv <- as.numeric(d[,var])
  
  # means by citizenship status
  m1 <- mean(d$dv[d$cit_citizenUS_alt == "naturalized US citizen"], na.rm=T)
  m2 <- mean(d$dv[d$cit_citizenUS_alt == "lawful permanent resident"], na.rm=T)
  m3 <- mean(d$dv[d$cit_citizenUS_alt == "likely unauthorized immigrant"], na.rm=T)
  
  # combine
  out <- data.frame(cbind(varname, m1,m2,m3))
  colnames(out) <- c("variable", "Naturalized", "LPR", "Unauth")
  out$Naturalized <- as.numeric(out$Naturalized)
  out$LPR <- as.numeric(out$LPR)
  out$Unauth <- as.numeric(out$Unauth)
  
  # output
  return(out)
}

### get distribution of immigration status first
aux1 <- data.frame(cbind("N", t(table(d$cit_citizenUS_alt[d$cit_citizenUS_alt != "natural-born US citizen"])[2:4])))
colnames(aux1) <- c("variable", "Naturalized", "LPR", "Unauth")
aux1$Naturalized <- as.numeric(aux1$Naturalized)
aux1$LPR <- as.numeric(aux1$LPR)
aux1$Unauth <- as.numeric(aux1$Unauth)

aux2 <- data.frame(cbind("Share", t(prop.table(table(d$cit_citizenUS_alt[d$cit_citizenUS_alt != "natural-born US citizen"])[2:4]))))
colnames(aux2) <- c("variable", "Naturalized", "LPR", "Unauth")
aux2$Naturalized <- as.numeric(aux2$Naturalized)
aux2$LPR <- as.numeric(aux2$LPR)
aux2$Unauth <- as.numeric(aux2$Unauth)


### get values
out <- rbind(aux1,aux2,
             cit.fun(var = "cit_citizenMexico", varname = "Mexican citizen (%)"),
             cit.fun(var = "language_english", varname = "English survey (%)"),
             cit.fun(var = "language_spanish", varname = "Spanish survey (%)"),
             cit.fun(var = "cit_yearsUS", varname = "Time in US (years)"),
             cit.fun(var = "bl_agree_immistatus_recode", varname = "My immigration status makes my life more difficult"),
             cit.fun(var = "bl_agree_deportation_recode", varname = "Even legal immigrants should be worried about deportation"),
             cit.fun(var = "bl_agree_leaveUS_recode", varname = "I sometimes think about leaving the US for good"),
             cit.fun(var = "bl_trust_ICE_dich", varname = "Trust in ICE (dichotomous)"),
             cit.fun(var = "cit_worries_deportation_dich", varname = "Worried about deportation"),
             cit.fun(var = "cit_spouse_citizenUS_citizen", varname = "Spouse is a US citizen"),
             cit.fun(var = "cit_spouse_citizenUS_unauth", varname = "Spouse is unauthorized"))

### Table output 
print(xtable(out, caption = "Comparison of select statistics by immigration status", align = c("lrccc")),
      caption.placement = "top", include.rownames=FALSE)





# Table A6: Randomization -------------------------------------------------

### Create output container
out <- data.frame(matrix(nrow = 14, ncol=3))
colnames(out) <- c("Experiment", "Attribute", "Share")
out$Experiment <- factor(c(1,1,2,2,3,3,3,3,4,4,4,4,5,5),
                         levels = 1:5,
                         labels = c("Overall", "Push", "Push - Positive", "Push - Negative", "Pull - additional questions"))
out$Attribute <- c("Push first", "Pull first",
                   "Positive change first", "Negative change first",
                   "Deportation", "Economic opportunities", "Healthcare", "Social network",
                   "Deportation", "Economic opportunities", "Healthcare", "Social network",
                   "Version 1", "Version 2")

### Get values
out[1:2,3] <- prop.table(table(d$first))
out[3:4,3] <- prop.table(table(d$pushSTART))
out[5:8,3] <- prop.table(table(d$pushPOSITIVE))
out[9:12,3] <- prop.table(table(d$pushNEGATIVE))
out[13:14,3] <- prop.table(table(d$E2_comparison))

### Table output 
print(xtable(out, caption = "Randomization", align = c("lrrc")),
      caption.placement = "top", include.rownames=FALSE)





# Figure A3: Attributes by order in which they were presented -------------

### create output container
out <- data.frame(matrix(nrow=49, ncol=4))
colnames(out) <- c("order", "dimension", "n", "share")
out$order <- rep(1:7, each=7)
out$dimension <- rep(levels(d$E2_dimension1), 7)

### get values
out$n[1:7] <- table(d$E2_dimension1)
out$n[8:14] <- table(d$E2_dimension2)
out$n[15:21] <- table(d$E2_dimension3)
out$n[22:28] <- table(d$E2_dimension4)
out$n[29:35] <- table(d$E2_dimension5)
out$n[36:42] <- table(d$E2_dimension6)
out$n[43:49] <- table(d$E2_dimension7)
out$share <- out$n / dim(d)[1]

### plot
p <- ggplot() +
  geom_bar(data = out, aes(x=share*100, y=dimension),
           stat="identity") + 
  xlab("Share of respondents (%)") +
  ylab("") +
  facet_wrap(~order) +
  theme_bw() + 
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title = element_text(size=14, face="bold"),
        axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=12),
        legend.position = "none")    
p
ggsave("FigureA3_Conjoint_Order.png", p, width=12,height=6)




# Figure A4: Distribution of levels ---------------------------------------

### output container
out <- data.frame(matrix(nrow=28, ncol=4))
colnames(out) <- c("dimension", "attribute", "n", "share")
out$dimension <-rep(levels(d$E2_dimension1), each=4)

### get values
out[1:4,2:4] <- cbind(levels(d.conjoint$E2_jobs), table(d.conjoint$E2_jobs),
                      prop.table(table(d.conjoint$E2_jobs)))
out[5:8,2:4] <- cbind(levels(d.conjoint$E2_rent), table(d.conjoint$E2_rent),
                      prop.table(table(d.conjoint$E2_rent)))
out[9:12,2:4] <- cbind(levels(d.conjoint$E2_netw), table(d.conjoint$E2_netw),
                       prop.table(table(d.conjoint$E2_netw)))
out[13:16,2:4] <- cbind(levels(d.conjoint$E2_immi), table(d.conjoint$E2_immi),
                        prop.table(table(d.conjoint$E2_immi)))
out[17:20,2:4] <- cbind(levels(d.conjoint$E2_deport), table(d.conjoint$E2_deport),
                        prop.table(table(d.conjoint$E2_deport)))
out[21:24,2:4] <- cbind(levels(d.conjoint$E2_health), table(d.conjoint$E2_health),
                        prop.table(table(d.conjoint$E2_health)))
out[25:28,2:4] <- cbind(levels(d.conjoint$E2_type), table(d.conjoint$E2_type),
                        prop.table(table(d.conjoint$E2_type)))
out$share <- as.numeric(out$share)

### plot
p <- ggplot() +
  geom_bar(data = out, aes(x=share*100, y=attribute),
           stat="identity") + 
  xlab("Share of respondents (%)") +
  scale_x_continuous(limits = c(0,30)) +
  ylab("") +
  facet_wrap(~dimension, scales="free_y") +
  theme_bw() + 
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title = element_text(size=14, face="bold"),
        axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=12),
        legend.position = "none")    
p
ggsave("FigureA4_Conjoint_Attributes.png", p, width=18,height=6)  






# Figure A5: Immigration status by treatment condition in the push --------

### Create output container
out <- data.frame(matrix(nrow=16, ncol=5))
colnames(out) <- c("class", "outcome", "treatment", "b", "se")
out$class <- factor(rep(rep(1:2,each=2),4), levels = 1:2,
                    labels = c("Deportation risk", "Healthcare access"))
out$outcome <- factor(rep(1:4,each=4), levels = 1:4,
                      labels = c("naturalized US citizen", "lawful permanent resident", "likely unauthorized immigrant", "complete response"))
out$treatment <- factor(rep(1:4,4), levels = 1:4,
                        labels = c("Higher\ndeportation risk", "Lower\ndeportation risk", 
                                   "Worse\nhealthcare access", "Better\nhealthcare access"))

### Get values
# naturalized US citizen
out[1,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="naturalized US citizen") ~ I(E1_treat == "policy: exclusionary (deport)") - 1, data=d.push))$coef[2,1:2]
out[2,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="naturalized US citizen") ~ I(E1_treat == "policy: welcoming (deport)") - 1, data=d.push))$coef[2,1:2]
out[3,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="naturalized US citizen") ~ I(E1_treat == "policy: exclusionary (health)") - 1, data=d.push))$coef[2,1:2]
out[4,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="naturalized US citizen") ~ I(E1_treat == "policy: welcoming (health)") - 1, data=d.push))$coef[2,1:2]

# lawful permanent resident
out[5,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="lawful permanent resident") ~ I(E1_treat == "policy: exclusionary (deport)") - 1, data=d.push))$coef[2,1:2]
out[6,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="lawful permanent resident") ~ I(E1_treat == "policy: welcoming (deport)") - 1, data=d.push))$coef[2,1:2]
out[7,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="lawful permanent resident") ~ I(E1_treat == "policy: exclusionary (health)") - 1, data=d.push))$coef[2,1:2]
out[8,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="lawful permanent resident") ~ I(E1_treat == "policy: welcoming (health)") - 1, data=d.push))$coef[2,1:2]

# likely unauthorized immigrant
out[9,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="likely unauthorized immigrant") ~ I(E1_treat == "policy: exclusionary (deport)") - 1, data=d.push))$coef[2,1:2]
out[10,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="likely unauthorized immigrant") ~ I(E1_treat == "policy: welcoming (deport)") - 1, data=d.push))$coef[2,1:2]
out[11,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="likely unauthorized immigrant") ~ I(E1_treat == "policy: exclusionary (health)") - 1, data=d.push))$coef[2,1:2]
out[12,4:5] <- summary(lm_robust(I(cit_citizenUS_alt=="likely unauthorized immigrant") ~ I(E1_treat == "policy: welcoming (health)") - 1, data=d.push))$coef[2,1:2]

# complete response
out[13,4:5] <- summary(lm_robust(I(complete=="complete response") ~ I(E1_treat == "policy: exclusionary (deport)") - 1, data=d.full.push))$coef[2,1:2]
out[14,4:5] <- summary(lm_robust(I(complete=="complete response") ~ I(E1_treat == "policy: welcoming (deport)") - 1, data=d.full.push))$coef[2,1:2]
out[15,4:5] <- summary(lm_robust(I(complete=="complete response") ~ I(E1_treat == "policy: exclusionary (health)") - 1, data=d.full.push))$coef[2,1:2]
out[16,4:5] <- summary(lm_robust(I(complete=="complete response") ~ I(E1_treat == "policy: welcoming (health)") - 1, data=d.full.push))$coef[2,1:2]


### Plot
p <- ggplot() + 
  geom_bar(data=out, aes(x=treatment, y=b*100, fill=class), stat="identity") +
  geom_errorbar(data=out, aes(x=treatment, ymin=(b-1.96*se)*100, ymax=(b+1.96*se)*100), 
                width=0.25, linewidth=1, color="gray30") +
  scale_fill_manual(values = c("burlywood3","burlywood4"), name="") +
  facet_wrap(~outcome) +
  xlab("Treatment group") +
  ylab("Share of respondents") +
  coord_flip() +
  theme_bw() +
  theme(strip.text = element_text(size=14,  face="bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size=12))   
p
ggsave("FigureA5_Push_impact_citizenship.png", p, width=8, height=6)




# Table A7: Summary statistics for moderators -----------------------------

### write function
summ.stats <- function(dv){
  d$dv <- as.numeric(d[,dv])
  mean <- mean(d$dv, na.rm=T)
  sd <- sd(d$dv, na.rm=T)
  return(cbind(mean,sd))
}

### output container
out <- data.frame(matrix(nrow=7, ncol=3))
colnames(out) <- c("Variable", "Mean", "SD")
out$Variable <- c("Has unauthorized family member",
                  "Knows someone who was deported",
                  "Trust in ICE officers",
                  "Worried about deportation",
                  "Immigration status makes life more difficult",
                  "Even legal immigrants should worry about deportation",
                  "Propensity to move away at baseline")

### get values
out[1,2:3] <- summ.stats(dv = "cit_family_unauthorized")
out[2,2:3] <- summ.stats(dv = "cit_know_deported")
out[3,2:3] <- summ.stats(dv = "bl_trust_ICE")
out[4,2:3] <- summ.stats(dv = "cit_worries_deportation")
out[5,2:3] <- summ.stats(dv = "bl_agree_immistatus")
out[6,2:3] <- summ.stats(dv = "bl_agree_deportation")
out[7,2:3] <- summ.stats(dv = "E1_dv.bl")


### output
print(xtable(out, caption = "Summary statistics", align = c("lrcc")),
      caption.placement = "top", include.rownames=FALSE)





# Table A8: Results of the push experiment --------------------------------

### Regressions
# baseline model
m1 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.full.push)

# baseline model -- complete responses
m2 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push)

# only naturalized citizens
m3 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt == "naturalized US citizen",])

# only LPRs
m4 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt == "lawful permanent resident",])

# only likely unauthorized immigrant
m5 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt == "likely unauthorized immigrant",])

# interaction
m6 <- felm(E1_dv ~ E1_treat*cit_citizenUS_alt | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen",])


### stargazer output
stargazer(m1,m2,m3,m4,m5,m6,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Pr(move away)"),
          column.labels = c("All", "All (complete)", "Naturalized", "LPR", "Unauthorized", "All (complete)"),
          covariate.labels=c("More economic opportunities",
                             "More personal connections", 
                             "Lower risk of deportation",
                             "Better access to healthcare",
                             "Fewer economic opportunities",
                             "Fewer personal connections",
                             "Higher risk of deportation",
                             "Worse access to healthcare",
                             "LPR",
                             "Likely unauthorized",
                             "LPR X More economic opportunities",
                             "LPR X More personal connections", 
                             "LPR X Lower risk of deportation",
                             "LPR X Better access to healthcare",
                             "LPR X Fewer economic opportunities",
                             "LPR X Fewer personal connections",
                             "LPR X Higher risk of deportation",
                             "LPR X Worse access to healthcare",
                             "Likely unauthorized X More economic opportunities",
                             "Likely unauthorized X More personal connections", 
                             "Likely unauthorized X Lower risk of deportation",
                             "Likely unauthorized X Better access to healthcare",
                             "Likely unauthorized X Fewer economic opportunities",
                             "Likely unauthorized X Fewer personal connections",
                             "Likely unauthorized X Higher risk of deportation",
                             "Likely unauthorized X Worse access to healthcare"),
          add.lines = list(
            c("Respondent-FE?", rep("YES", 6))), 
          digits = 3)




# Table A9: Influence of select moderators --------------------------------

### Regressions
d.push$moderator <- d.push$bl_agree_deportation_dich
m1 <- felm(E1_dv ~ E1_treat*moderator | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")
m2 <- felm(E1_dv ~ E1_treat*moderator + E1_treat*noncit_dich | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")

d.push$moderator <- d.push$bl_agree_immistatus_dich
m3 <- felm(E1_dv ~ E1_treat*moderator | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")
m4 <- felm(E1_dv ~ E1_treat*moderator + E1_treat*noncit_dich | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")

d.push$moderator <- d.push$cit_family_unauthorized
m5 <- felm(E1_dv ~ E1_treat*moderator | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")
m6 <- felm(E1_dv ~ E1_treat*moderator + E1_treat*noncit_dich | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")

d.push$moderator <- d.push$cit_know_deported
m7 <- felm(E1_dv ~ E1_treat*moderator | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")
m8 <- felm(E1_dv ~ E1_treat*moderator + E1_treat*noncit_dich | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen")


### stargazer output
stargazer(m1,m2,m3,m4,m5,m6,m7,m8,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Pr(move away)"),
          column.labels = c("Worried about deportation", "Even legal immigrants", "Unauth. family", "Know deportee"),
          column.separate = c(2, 2, 2, 2),
          covariate.labels=c("More economic opportunities",
                             "More personal connections", 
                             "Lower risk of deportation",
                             "Better access to healthcare",
                             "Fewer economic opportunities",
                             "Fewer personal connections",
                             "Higher risk of deportation",
                             "Worse access to healthcare",
                             "Moderator",
                             "Non-citizen",
                             "Moderator X More economic opportunities",
                             "Moderator X More personal connections", 
                             "Moderator X Lower risk of deportation",
                             "Moderator X Better access to healthcare",
                             "Moderator X Fewer economic opportunities",
                             "Moderator X Fewer personal connections",
                             "Moderator X Higher risk of deportation",
                             "Moderator X Worse access to healthcare",
                             "Non-citizen X More economic opportunities",
                             "Non-citizen X More personal connections", 
                             "Non-citizen X Lower risk of deportation",
                             "Non-citizen X Better access to healthcare",
                             "Non-citizen X Fewer economic opportunities",
                             "Non-citizen X Fewer personal connections",
                             "Non-citizen X Higher risk of deportation",
                             "Non-citizen X Worse access to healthcare"),
          add.lines = list(
            c("Respondent-FE?", rep("YES", 8))), 
          digits = 3)





# Table A10: Interaction with state immigration environment ---------------

### Regressions
# Democratic governor
d.push$moderator <- I(d.push$DemGovernor)
m1 <- felm(E1_dv ~ E1_treat*moderator | ResponseId | 0 | ResponseId, data=d.push)

d.push$moderator <- d.push$DemLegislature
m2 <- felm(E1_dv ~ E1_treat*moderator | ResponseId | 0 | ResponseId, data=d.push)


### stargazer output
stargazer(m1,m2,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Pr(move away)"),
          #column.labels = c("Welcoming environment (Monogan)", "Welcoming environment (Urban Institute)", "Democratic Governor", "Democratic State Legislature"),
          column.labels = c("Democratic Governor", "Democratic State Legislature"),
          covariate.labels=c("More economic opportunities",
                             "More personal connections", 
                             "Lower risk of deportation",
                             "Better access to healthcare",
                             "Fewer economic opportunities",
                             "Fewer personal connections",
                             "Higher risk of deportation",
                             "Worse access to healthcare",
                             "Moderator",
                             "Moderator X More economic opportunities",
                             "Moderator X More personal connections", 
                             "Moderator X Lower risk of deportation",
                             "Moderator X Better access to healthcare",
                             "Moderator X Fewer economic opportunities",
                             "Moderator X Fewer personal connections",
                             "Moderator X Higher risk of deportation",
                             "Moderator X Worse access to healthcare"),
          add.lines = list(
            c("Respondent-FE?", rep("YES", 2))), 
          digits = 3)




# Figure A11: Summary of robustness checks --------------------------------

### Regressions
# only push first
m1 <- felm(E1_dv ~ E1_treat*cit_citizenUS_alt | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen" & 
             d.push$first == "push first")

# only first treatment group
m2 <- felm(E1_dv ~ E1_treat*cit_citizenUS_alt | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen" & 
             d.push$E1_period != 3)

# only Mexican citizens
m3 <- felm(E1_dv ~ E1_treat*cit_citizenUS_alt | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen" & 
             d.push$cit_citizenMexico == 1)

# passed attention check
m4 <- felm(E1_dv ~ E1_treat*cit_citizenUS_alt | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen" & 
             d.push$attentioncheck == 1)

# duration
m5 <- felm(E1_dv ~ E1_treat*cit_citizenUS_alt | ResponseId | 0 | ResponseId, data=d.push,
           subset=d.push$cit_citizenUS_alt != "natural-born US citizen" & 
             d.push$duration >= 0.5*median(d.push$duration, na.rm=T) & d.push$duration <= 2*median(d.push$duration, na.rm=T))




### stargazer output
stargazer(m1,m2,m3,m4,m5,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Pr(move away)"),
          column.labels = c("Push first", "First treatment", "Mexican citizen", "Attention check", "Duration"),
          covariate.labels=c("More economic opportunities",
                             "More personal connections", 
                             "Lower risk of deportation",
                             "Better access to healthcare",
                             "Fewer economic opportunities",
                             "Fewer personal connections",
                             "Higher risk of deportation",
                             "Worse access to healthcare",
                             "LPR",
                             "Likely unauthorized",
                             "LPR X More economic opportunities",
                             "LPR X More personal connections", 
                             "LPR X Lower risk of deportation",
                             "LPR X Better access to healthcare",
                             "LPR X Fewer economic opportunities",
                             "LPR X Fewer personal connections",
                             "LPR X Higher risk of deportation",
                             "LPR X Worse access to healthcare",
                             "Likely unauthorized X More economic opportunities",
                             "Likely unauthorized X More personal connections", 
                             "Likely unauthorized X Lower risk of deportation",
                             "Likely unauthorized X Better access to healthcare",
                             "Likely unauthorized X Fewer economic opportunities",
                             "Likely unauthorized X Fewer personal connections",
                             "Likely unauthorized X Higher risk of deportation",
                             "Likely unauthorized X Worse access to healthcare"),
          add.lines = list(
            c("Respondent-FE?", rep("YES", 5))), 
          digits = 3)




# Figure A6: Differences in marginal means 1: moderators ------------------


### Have unauthorized family member
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~cit_family_unauthorized)

# prepare output container
out1 <- data.frame(matrix(nrow=35, ncol=6))
colnames(out1) <- c("moderator","dimension", "attribute", "b", "lb", "ub")
out1$moderator <- 3 #"Unauthorized\nfamily member"
out1$dimension <- rep(1:7, each=5)
out1$dimension <- factor(out1$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out1$attribute <- 35:1
out1$attribute <- factor(out1$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out1[2:5,4:6] <- reg[1:4,c(6,10,11)]
out1[7:10,4:6] <- reg[5:8,c(6,10,11)]
out1[12:15,4:6] <- reg[9:12,c(6,10,11)]
out1[17:20,4:6] <- reg[13:16,c(6,10,11)]
out1[22:25,4:6] <- reg[17:20,c(6,10,11)]
out1[27:30,4:6] <- reg[21:24,c(6,10,11)]
out1[32:35,4:6] <- reg[25:28,c(6,10,11)]



### Know someone who was deported
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~cit_know_deported)

# prepare output container
out2 <- data.frame(matrix(nrow=35, ncol=6))
colnames(out2) <- c("moderator","dimension", "attribute", "b", "lb", "ub")
out2$moderator <- 4 #"Knows someone\nwho was deported"
out2$dimension <- rep(1:7, each=5)
out2$dimension <- factor(out2$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out2$attribute <- 35:1
out2$attribute <- factor(out2$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out2[2:5,4:6] <- reg[1:4,c(6,10,11)]
out2[7:10,4:6] <- reg[5:8,c(6,10,11)]
out2[12:15,4:6] <- reg[9:12,c(6,10,11)]
out2[17:20,4:6] <- reg[13:16,c(6,10,11)]
out2[22:25,4:6] <- reg[17:20,c(6,10,11)]
out2[27:30,4:6] <- reg[21:24,c(6,10,11)]
out2[32:35,4:6] <- reg[25:28,c(6,10,11)]


### Even legal immigrants should be worried about deportation
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~bl_agree_deportation_dich)

# prepare output container
out3 <- data.frame(matrix(nrow=35, ncol=6))
colnames(out3) <- c("moderator","dimension", "attribute", "b", "lb", "ub")
out3$moderator <- 2 #"Even legal immigrants\nshould worry about deportation"
out3$dimension <- rep(1:7, each=5)
out3$dimension <- factor(out3$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out3$attribute <- 35:1
out3$attribute <- factor(out3$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out3[2:5,4:6] <- reg[1:4,c(6,10,11)]
out3[7:10,4:6] <- reg[5:8,c(6,10,11)]
out3[12:15,4:6] <- reg[9:12,c(6,10,11)]
out3[17:20,4:6] <- reg[13:16,c(6,10,11)]
out3[22:25,4:6] <- reg[17:20,c(6,10,11)]
out3[27:30,4:6] <- reg[21:24,c(6,10,11)]
out3[32:35,4:6] <- reg[25:28,c(6,10,11)]


### Worried about deportation
# regression
reg <- cj(d.conjoint, 
          E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~cit_worries_deportation_dich)

# prepare output container
out4 <- data.frame(matrix(nrow=35, ncol=6))
colnames(out4) <- c("moderator","dimension", "attribute", "b", "lb", "ub")
out4$moderator <- 1 #"Even legal immigrants\nshould worry about deportation"
out4$dimension <- rep(1:7, each=5)
out4$dimension <- factor(out4$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out4$attribute <- 35:1
out4$attribute <- factor(out4$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out4[2:5,4:6] <- reg[1:4,c(6,10,11)]
out4[7:10,4:6] <- reg[5:8,c(6,10,11)]
out4[12:15,4:6] <- reg[9:12,c(6,10,11)]
out4[17:20,4:6] <- reg[13:16,c(6,10,11)]
out4[22:25,4:6] <- reg[17:20,c(6,10,11)]
out4[27:30,4:6] <- reg[21:24,c(6,10,11)]
out4[32:35,4:6] <- reg[25:28,c(6,10,11)]



### combine all four
out <- rbind(out1,out2)
out <- rbind(out,out3)
out <- rbind(out,out4)
out$sample <- "Yes - No"


### factorize moderators
out$moderator <- factor(out$moderator, levels = 1:4,
                        labels = c("Worried\nabout deportation",
                                   "Even legal immigrants\nshould worry\nabout deportation",
                                   "Unauthorized\nfamily member",
                                   "Knows someone\nwho was deported"))


### plot
p <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$dimension %in% c("Deportation risk"),], 
             aes(x=attribute, y=b*100, col=sample), size=2.5,
             position = position_dodge(width=0.5)) +
  geom_errorbar(data=out[out$dimension %in% c("Deportation risk"),], 
                aes(x=attribute, ymin=lb*100, ymax=ub*100, col=sample), width=0,
                position = position_dodge(width=0.5),size=.75) + 
  facet_wrap(~moderator, nrow=1, scales="free_x") +
  scale_color_manual(values = c("burlywood4")) +
  coord_flip() + 
  xlab("") +
  ylab("Difference in Marginal Means") + 
  theme_bw() +  
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_text(size=14, face="bold"),
        axis.text.y = element_text(face=rep(c(rep("plain",4), "bold"),2), size=12),
        axis.text.x = element_text(size=12),
        axis.ticks.y = element_blank(),
        strip.text = element_text(size=14, face="bold"))
p
ggsave("FigureA6_Pull_results_HTE_MMDiff.png", p, width=13, height=4)





# Figure A7: Differences in marginal means 2: other outcomes --------------

### Fear discrimination
# regression
reg <- cj(d.conjoint, 
          E2_discrimination ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~cit_citizenUS_altDICH)

# prepare output container
out1 <- data.frame(matrix(nrow=35, ncol=7))
colnames(out1) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out1$outcome <- "Fear\ndiscrimination"
out1$sample <- 1
out1$dimension <- rep(1:7, each=5)
out1$dimension <- factor(out1$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out1$attribute <- 35:1
out1$attribute <- factor(out1$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out1[2:5,5:7] <- reg[1:4,c(6,10,11)]
out1[7:10,5:7] <- reg[5:8,c(6,10,11)]
out1[12:15,5:7] <- reg[9:12,c(6,10,11)]
out1[17:20,5:7] <- reg[13:16,c(6,10,11)]
out1[22:25,5:7] <- reg[17:20,c(6,10,11)]
out1[27:30,5:7] <- reg[21:24,c(6,10,11)]
out1[32:35,5:7] <- reg[25:28,c(6,10,11)]



### Be safe
# regression
reg <- cj(d.conjoint, 
          E2_besafe ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~cit_citizenUS_altDICH)

# prepare output container
out2 <- data.frame(matrix(nrow=35, ncol=7))
colnames(out2) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out2$outcome <- "Be safe"
out2$sample <- 1
out2$dimension <- rep(1:7, each=5)
out2$dimension <- factor(out2$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out2$attribute <- 35:1
out2$attribute <- factor(out2$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out2[2:5,5:7] <- reg[1:4,c(6,10,11)]
out2[7:10,5:7] <- reg[5:8,c(6,10,11)]
out2[12:15,5:7] <- reg[9:12,c(6,10,11)]
out2[17:20,5:7] <- reg[13:16,c(6,10,11)]
out2[22:25,5:7] <- reg[17:20,c(6,10,11)]
out2[27:30,5:7] <- reg[21:24,c(6,10,11)]
out2[32:35,5:7] <- reg[25:28,c(6,10,11)]



### Local government protects immigrants
# regression
reg <- cj(d.conjoint, 
          E2_govprotect ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~cit_citizenUS_altDICH)

# prepare output container
out3 <- data.frame(matrix(nrow=35, ncol=7))
colnames(out3) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out3$outcome <- "Gov't protects\nimmigrants"
out3$sample <- 1
out3$dimension <- rep(1:7, each=5)
out3$dimension <- factor(out3$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out3$attribute <- 35:1
out3$attribute <- factor(out3$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out3[2:5,5:7] <- reg[1:4,c(6,10,11)]
out3[7:10,5:7] <- reg[5:8,c(6,10,11)]
out3[12:15,5:7] <- reg[9:12,c(6,10,11)]
out3[17:20,5:7] <- reg[13:16,c(6,10,11)]
out3[22:25,5:7] <- reg[17:20,c(6,10,11)]
out3[27:30,5:7] <- reg[21:24,c(6,10,11)]
out3[32:35,5:7] <- reg[25:28,c(6,10,11)]


### Government is liberal
# regression
reg <- cj(d.conjoint, 
          E2_govliberal ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_difference",
          by=~cit_citizenUS_altDICH)

# prepare output container
out4 <- data.frame(matrix(nrow=35, ncol=7))
colnames(out4) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out4$outcome <- "Gov't is\nliberal"
out4$sample <- 1
out4$dimension <- rep(1:7, each=5)
out4$dimension <- factor(out4$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out4$attribute <- 35:1
out4$attribute <- factor(out4$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out4[2:5,5:7] <- reg[1:4,c(6,10,11)]
out4[7:10,5:7] <- reg[5:8,c(6,10,11)]
out4[12:15,5:7] <- reg[9:12,c(6,10,11)]
out4[17:20,5:7] <- reg[13:16,c(6,10,11)]
out4[22:25,5:7] <- reg[17:20,c(6,10,11)]
out4[27:30,5:7] <- reg[21:24,c(6,10,11)]
out4[32:35,5:7] <- reg[25:28,c(6,10,11)]



### combine all four
out <- rbind(out1,out2)
out <- rbind(out,out3)
out <- rbind(out,out4)

out$sample <- "Likely unauth. - naturalized"


### Plot
p <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out[out$dimension %in% c("Deportation risk"),], 
             aes(x=attribute, y=b*100, shape=sample, col=sample), size=2.5,
             position = position_dodge(width=0.5)) +
  geom_errorbar(data=out[out$dimension %in% c("Deportation risk"),], 
                aes(x=attribute, ymin=lb*100, ymax=ub*100, col=sample), width=0,
                position = position_dodge(width=0.5),size=.75) + 
  scale_color_manual(values = c("burlywood4"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  scale_linetype_manual(values=c(2,1), name="") +
  #scale_y_continuous(breaks = seq(20,80,10)) +
  facet_wrap(~outcome, nrow=1, scales="free_x") +
  coord_flip() + 
  xlab("") +
  ylab("Difference in Marginal Means") + 
  theme_bw() +  
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_text(size=14, face="bold"),
        axis.text.y = element_text(face=rep(c(rep("plain",4), "bold"),2), size=12),
        axis.text.x = element_text(size=12),
        axis.ticks.y = element_blank(),
        strip.text = element_text(size=14, face="bold"))
p
ggsave("FigureA7_Pull_results_other_MMDiff.png", p, width=11, height=3.5)





# Figure A8: Conjoint results when estimating AMCEs -----------------------

### Plot 1: all respondents
# regression
reg <- amce(data = d.conjoint, formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- reg[1:4,c(5,9,10)]
out[7:10,3:5] <- reg[5:8,c(5,9,10)]
out[12:15,3:5] <- reg[9:12,c(5,9,10)]
out[17:20,3:5] <- reg[13:16,c(5,9,10)]
out[22:25,3:5] <- reg[17:20,c(5,9,10)]
out[27:30,3:5] <- reg[21:24,c(5,9,10)]
out[32:35,3:5] <- reg[25:28,c(5,9,10)]

# define sample
out$sample <- "all respondents"

# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2) +
  #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0) + 
  scale_color_manual(values = "chocolate4", name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("AMCE") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


### Plot 2: by citizenship
# Regression
mm_by <- cj(d.conjoint, 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "amce",
            by=~cit_citizenUS_alt)

# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("AMCE") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


### combining all three plots
p <- ggarrange(p1,p2,widths=c(0.45,0.25), nrow = 1)
p
ggsave("FigureA8_Pull_results_AMCE.png", p, width=10, height=7)




# Figure A9: Secondary outcomes -------------------------------------------

### Be economically secure
# regression
reg <- cj(d.conjoint, 
          E2_econsecure ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out1 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out1) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out1$outcome <- "Be economically\nsecure"
out1$sample <- rep(1:2,each=35)
out1$sample <- factor(out1$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out1$dimension <- rep(rep(1:7, each=5),2)
out1$dimension <- factor(out1$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out1$attribute <- rep(35:1,2)
out1$attribute <- factor(out1$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out1[2:5,5:7] <- reg[1:4,c(6,10,11)]
out1[7:10,5:7] <- reg[5:8,c(6,10,11)]
out1[12:15,5:7] <- reg[9:12,c(6,10,11)]
out1[17:20,5:7] <- reg[13:16,c(6,10,11)]
out1[22:25,5:7] <- reg[17:20,c(6,10,11)]
out1[27:30,5:7] <- reg[21:24,c(6,10,11)]
out1[32:35,5:7] <- reg[25:28,c(6,10,11)]

out1[37:40,5:7] <- reg[29:32,c(6,10,11)]
out1[42:45,5:7] <- reg[33:36,c(6,10,11)]
out1[47:50,5:7] <- reg[37:40,c(6,10,11)]
out1[52:55,5:7] <- reg[41:44,c(6,10,11)]
out1[57:60,5:7] <- reg[45:48,c(6,10,11)]
out1[62:65,5:7] <- reg[49:52,c(6,10,11)]
out1[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Have many opportunities
# regression
reg <- cj(d.conjoint, 
          E2_opportunities ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out2 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out2) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out2$outcome <- "Have many\nopportunities"
out2$sample <- rep(1:2,each=35)
out2$sample <- factor(out2$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out2$dimension <- rep(rep(1:7, each=5),2)
out2$dimension <- factor(out2$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out2$attribute <- rep(35:1,2)
out2$attribute <- factor(out2$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out2[2:5,5:7] <- reg[1:4,c(6,10,11)]
out2[7:10,5:7] <- reg[5:8,c(6,10,11)]
out2[12:15,5:7] <- reg[9:12,c(6,10,11)]
out2[17:20,5:7] <- reg[13:16,c(6,10,11)]
out2[22:25,5:7] <- reg[17:20,c(6,10,11)]
out2[27:30,5:7] <- reg[21:24,c(6,10,11)]
out2[32:35,5:7] <- reg[25:28,c(6,10,11)]

out2[37:40,5:7] <- reg[29:32,c(6,10,11)]
out2[42:45,5:7] <- reg[33:36,c(6,10,11)]
out2[47:50,5:7] <- reg[37:40,c(6,10,11)]
out2[52:55,5:7] <- reg[41:44,c(6,10,11)]
out2[57:60,5:7] <- reg[45:48,c(6,10,11)]
out2[62:65,5:7] <- reg[49:52,c(6,10,11)]
out2[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Have many people to help
# regression
reg <- cj(d.conjoint, 
          E2_help ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out3 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out3) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out3$outcome <- "Have many people\nwho can help"
out3$sample <- rep(1:2,each=35)
out3$sample <- factor(out3$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out3$dimension <- rep(rep(1:7, each=5),2)
out3$dimension <- factor(out3$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out3$attribute <- rep(35:1,2)
out3$attribute <- factor(out3$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out3[2:5,5:7] <- reg[1:4,c(6,10,11)]
out3[7:10,5:7] <- reg[5:8,c(6,10,11)]
out3[12:15,5:7] <- reg[9:12,c(6,10,11)]
out3[17:20,5:7] <- reg[13:16,c(6,10,11)]
out3[22:25,5:7] <- reg[17:20,c(6,10,11)]
out3[27:30,5:7] <- reg[21:24,c(6,10,11)]
out3[32:35,5:7] <- reg[25:28,c(6,10,11)]

out3[37:40,5:7] <- reg[29:32,c(6,10,11)]
out3[42:45,5:7] <- reg[33:36,c(6,10,11)]
out3[47:50,5:7] <- reg[37:40,c(6,10,11)]
out3[52:55,5:7] <- reg[41:44,c(6,10,11)]
out3[57:60,5:7] <- reg[45:48,c(6,10,11)]
out3[62:65,5:7] <- reg[49:52,c(6,10,11)]
out3[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Feel alone
# regression
reg <- cj(d.conjoint, 
          E2_alone ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out4 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out4) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out4$outcome <- "Feel alone"
out4$sample <- rep(1:2,each=35)
out4$sample <- factor(out4$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out4$dimension <- rep(rep(1:7, each=5),2)
out4$dimension <- factor(out4$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out4$attribute <- rep(35:1,2)
out4$attribute <- factor(out4$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out4[2:5,5:7] <- reg[1:4,c(6,10,11)]
out4[7:10,5:7] <- reg[5:8,c(6,10,11)]
out4[12:15,5:7] <- reg[9:12,c(6,10,11)]
out4[17:20,5:7] <- reg[13:16,c(6,10,11)]
out4[22:25,5:7] <- reg[17:20,c(6,10,11)]
out4[27:30,5:7] <- reg[21:24,c(6,10,11)]
out4[32:35,5:7] <- reg[25:28,c(6,10,11)]

out4[37:40,5:7] <- reg[29:32,c(6,10,11)]
out4[42:45,5:7] <- reg[33:36,c(6,10,11)]
out4[47:50,5:7] <- reg[37:40,c(6,10,11)]
out4[52:55,5:7] <- reg[41:44,c(6,10,11)]
out4[57:60,5:7] <- reg[45:48,c(6,10,11)]
out4[62:65,5:7] <- reg[49:52,c(6,10,11)]
out4[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Be safe
# regression
reg <- cj(d.conjoint, 
          E2_besafe ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out5 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out5) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out5$outcome <- "Be safe"
out5$sample <- rep(1:2,each=35)
out5$sample <- factor(out5$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out5$dimension <- rep(rep(1:7, each=5),2)
out5$dimension <- factor(out5$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out5$attribute <- rep(35:1,2)
out5$attribute <- factor(out5$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out5[2:5,5:7] <- reg[1:4,c(6,10,11)]
out5[7:10,5:7] <- reg[5:8,c(6,10,11)]
out5[12:15,5:7] <- reg[9:12,c(6,10,11)]
out5[17:20,5:7] <- reg[13:16,c(6,10,11)]
out5[22:25,5:7] <- reg[17:20,c(6,10,11)]
out5[27:30,5:7] <- reg[21:24,c(6,10,11)]
out5[32:35,5:7] <- reg[25:28,c(6,10,11)]

out5[37:40,5:7] <- reg[29:32,c(6,10,11)]
out5[42:45,5:7] <- reg[33:36,c(6,10,11)]
out5[47:50,5:7] <- reg[37:40,c(6,10,11)]
out5[52:55,5:7] <- reg[41:44,c(6,10,11)]
out5[57:60,5:7] <- reg[45:48,c(6,10,11)]
out5[62:65,5:7] <- reg[49:52,c(6,10,11)]
out5[67:70,5:7] <- reg[53:56,c(6,10,11)]



### Local government protects immigrants
# regression
reg <- cj(d.conjoint, 
          E2_govprotect ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out6 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out6) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out6$outcome <- "Gov't protects\nimmigrants"
out6$sample <- rep(1:2,each=35)
out6$sample <- factor(out6$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out6$dimension <- rep(rep(1:7, each=5),2)
out6$dimension <- factor(out6$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out6$attribute <- rep(35:1,2)
out6$attribute <- factor(out6$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out6[2:5,5:7] <- reg[1:4,c(6,10,11)]
out6[7:10,5:7] <- reg[5:8,c(6,10,11)]
out6[12:15,5:7] <- reg[9:12,c(6,10,11)]
out6[17:20,5:7] <- reg[13:16,c(6,10,11)]
out6[22:25,5:7] <- reg[17:20,c(6,10,11)]
out6[27:30,5:7] <- reg[21:24,c(6,10,11)]
out6[32:35,5:7] <- reg[25:28,c(6,10,11)]

out6[37:40,5:7] <- reg[29:32,c(6,10,11)]
out6[42:45,5:7] <- reg[33:36,c(6,10,11)]
out6[47:50,5:7] <- reg[37:40,c(6,10,11)]
out6[52:55,5:7] <- reg[41:44,c(6,10,11)]
out6[57:60,5:7] <- reg[45:48,c(6,10,11)]
out6[62:65,5:7] <- reg[49:52,c(6,10,11)]
out6[67:70,5:7] <- reg[53:56,c(6,10,11)]



### Government is liberal
# regression
reg <- cj(d.conjoint, 
          E2_govliberal ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out7 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out7) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out7$outcome <- "Gov't is\nliberal"
out7$sample <- rep(1:2,each=35)
out7$sample <- factor(out7$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out7$dimension <- rep(rep(1:7, each=5),2)
out7$dimension <- factor(out7$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out7$attribute <- rep(35:1,2)
out7$attribute <- factor(out7$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out7[2:5,5:7] <- reg[1:4,c(6,10,11)]
out7[7:10,5:7] <- reg[5:8,c(6,10,11)]
out7[12:15,5:7] <- reg[9:12,c(6,10,11)]
out7[17:20,5:7] <- reg[13:16,c(6,10,11)]
out7[22:25,5:7] <- reg[17:20,c(6,10,11)]
out7[27:30,5:7] <- reg[21:24,c(6,10,11)]
out7[32:35,5:7] <- reg[25:28,c(6,10,11)]

out7[37:40,5:7] <- reg[29:32,c(6,10,11)]
out7[42:45,5:7] <- reg[33:36,c(6,10,11)]
out7[47:50,5:7] <- reg[37:40,c(6,10,11)]
out7[52:55,5:7] <- reg[41:44,c(6,10,11)]
out7[57:60,5:7] <- reg[45:48,c(6,10,11)]
out7[62:65,5:7] <- reg[49:52,c(6,10,11)]
out7[67:70,5:7] <- reg[53:56,c(6,10,11)]



### Trust the police
# regression
reg <- cj(d.conjoint, 
          E2_trustpolice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out8 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out8) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out8$outcome <- "Trust the police"
out8$sample <- rep(1:2,each=35)
out8$sample <- factor(out8$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out8$dimension <- rep(rep(1:7, each=5),2)
out8$dimension <- factor(out8$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out8$attribute <- rep(35:1,2)
out8$attribute <- factor(out8$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out8[2:5,5:7] <- reg[1:4,c(6,10,11)]
out8[7:10,5:7] <- reg[5:8,c(6,10,11)]
out8[12:15,5:7] <- reg[9:12,c(6,10,11)]
out8[17:20,5:7] <- reg[13:16,c(6,10,11)]
out8[22:25,5:7] <- reg[17:20,c(6,10,11)]
out8[27:30,5:7] <- reg[21:24,c(6,10,11)]
out8[32:35,5:7] <- reg[25:28,c(6,10,11)]

out8[37:40,5:7] <- reg[29:32,c(6,10,11)]
out8[42:45,5:7] <- reg[33:36,c(6,10,11)]
out8[47:50,5:7] <- reg[37:40,c(6,10,11)]
out8[52:55,5:7] <- reg[41:44,c(6,10,11)]
out8[57:60,5:7] <- reg[45:48,c(6,10,11)]
out8[62:65,5:7] <- reg[49:52,c(6,10,11)]
out8[67:70,5:7] <- reg[53:56,c(6,10,11)]


### High quality of life
# regression
reg <- cj(d.conjoint, 
          E2_quallife ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out9 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out9) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out9$outcome <- "High quality\nof life"
out9$sample <- rep(1:2,each=35)
out9$sample <- factor(out9$sample, levels = 1:2,
                      labels = c("naturalized", "likely unauthorized"))
out9$dimension <- rep(rep(1:7, each=5),2)
out9$dimension <- factor(out9$dimension, levels=1:7, 
                         labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                    "Social network", "Share Mexican immigrants",
                                    "Deportation risk", "Healthcare access",
                                    "Type of community"))
out9$attribute <- rep(35:1,2)
out9$attribute <- factor(out9$attribute, levels=1:35,
                         labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                    "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                    "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                    "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                    "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                    "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                    "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                    "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out9[2:5,5:7] <- reg[1:4,c(6,10,11)]
out9[7:10,5:7] <- reg[5:8,c(6,10,11)]
out9[12:15,5:7] <- reg[9:12,c(6,10,11)]
out9[17:20,5:7] <- reg[13:16,c(6,10,11)]
out9[22:25,5:7] <- reg[17:20,c(6,10,11)]
out9[27:30,5:7] <- reg[21:24,c(6,10,11)]
out9[32:35,5:7] <- reg[25:28,c(6,10,11)]

out9[37:40,5:7] <- reg[29:32,c(6,10,11)]
out9[42:45,5:7] <- reg[33:36,c(6,10,11)]
out9[47:50,5:7] <- reg[37:40,c(6,10,11)]
out9[52:55,5:7] <- reg[41:44,c(6,10,11)]
out9[57:60,5:7] <- reg[45:48,c(6,10,11)]
out9[62:65,5:7] <- reg[49:52,c(6,10,11)]
out9[67:70,5:7] <- reg[53:56,c(6,10,11)]


### Fear discrimination
# regression
reg <- cj(d.conjoint, 
          E2_discrimination ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
          by=~cit_citizenUS_altDICH)

# prepare output container
out10 <- data.frame(matrix(nrow=70, ncol=7))
colnames(out10) <- c("outcome","sample","dimension", "attribute", "b", "lb", "ub")
out10$outcome <- "Fear\ndiscrimination"
out10$sample <- rep(1:2,each=35)
out10$sample <- factor(out10$sample, levels = 1:2,
                       labels = c("naturalized", "likely unauthorized"))
out10$dimension <- rep(rep(1:7, each=5),2)
out10$dimension <- factor(out10$dimension, levels=1:7, 
                          labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                     "Social network", "Share Mexican immigrants",
                                     "Deportation risk", "Healthcare access",
                                     "Type of community"))
out10$attribute <- rep(35:1,2)
out10$attribute <- factor(out10$attribute, levels=1:35,
                          labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                     "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                     "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                     "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                     "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                     "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                     "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                     "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out10[2:5,5:7] <- reg[1:4,c(6,10,11)]
out10[7:10,5:7] <- reg[5:8,c(6,10,11)]
out10[12:15,5:7] <- reg[9:12,c(6,10,11)]
out10[17:20,5:7] <- reg[13:16,c(6,10,11)]
out10[22:25,5:7] <- reg[17:20,c(6,10,11)]
out10[27:30,5:7] <- reg[21:24,c(6,10,11)]
out10[32:35,5:7] <- reg[25:28,c(6,10,11)]

out10[37:40,5:7] <- reg[29:32,c(6,10,11)]
out10[42:45,5:7] <- reg[33:36,c(6,10,11)]
out10[47:50,5:7] <- reg[37:40,c(6,10,11)]
out10[52:55,5:7] <- reg[41:44,c(6,10,11)]
out10[57:60,5:7] <- reg[45:48,c(6,10,11)]
out10[62:65,5:7] <- reg[49:52,c(6,10,11)]
out10[67:70,5:7] <- reg[53:56,c(6,10,11)]


### combine all ten
out <- rbind(out1,out2)
out <- rbind(out,out3)
out <- rbind(out,out4)
out <- rbind(out,out5)
out <- rbind(out,out6)
out <- rbind(out,out7)
out <- rbind(out,out8)
out <- rbind(out,out9)
out <- rbind(out,out10)


### Plot
p <- ggplot() + 
  geom_hline(yintercept = 50, lty=2) +
  geom_point(data=out, 
             aes(x=attribute, y=b*100, shape=sample, col=sample), size=2.5,
             position = position_dodge(width=0.5)) +
  geom_errorbar(data=out, 
                aes(x=attribute, ymin=lb*100, ymax=ub*100, col=sample), width=0,
                position = position_dodge(width=0.5),size=.75) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values=c(19,17), name="") +
  scale_y_continuous(breaks = seq(20,80,10)) +
  facet_wrap(~outcome, nrow=2, scales="free_x") +
  coord_flip() + 
  xlab("") +
  ylab("Marginal Mean") + 
  theme_bw() +  
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title = element_text(size=14, face="bold"),
        axis.text.y = element_text(face=rep(c(rep("plain",4), "bold"),2), size=12),
        axis.text.x = element_text(size=12),
        axis.ticks.y = element_blank(),
        strip.text = element_text(size=14, face="bold"))
p
ggsave("FigureA9_Pull_results_other_all.png", p, width=13, height=15)





# Figure A10: Conjoint results when restricting the analysis to re --------


### subset data
d.conjoint2 <- subset(d.conjoint, d.conjoint$first == "pull first")

### Regression
mm_by <- cj(d.conjoint2, 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

### prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

### define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

### get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

### plot
p <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p
ggsave("FigureA10_Pull_robustness_pullfirst.png", p, width=8, height=7)





# Figure A11: Conjoint results when restricting the analysis to th --------

### subset data
d.conjoint2 <- subset(d.conjoint, d.conjoint$E2_city == "A" | d.conjoint$E2_city == "B")

### Regression
mm_by <- cj(d.conjoint2, 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

### prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

### define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

### get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

### plot
p <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p
ggsave("FigureA11_Pull_robustness_firstpair.png", p, width=8, height=7)





# Figure A12: Conjoint results when restricting the analysis to Me --------

### subset data
d.conjoint2 <- subset(d.conjoint, d.conjoint$cit_citizenMexico == 1)

### Regression
mm_by <- cj(d.conjoint2, 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

### prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

### define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

### get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

### plot
p <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p
ggsave("FigureA12_Pull_robustness_Mexican.png", p, width=8, height=7)



# Figure A13: Conjoint results when dropping respondents who faile --------

### subset data
d.conjoint2 <- subset(d.conjoint, d.conjoint$attentioncheck == 1)

### Regression
mm_by <- cj(d.conjoint2, 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

### prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

### define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

### get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

### plot
p <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p
ggsave("FigureA13_Pull_robustness_attentioncheck.png", p, width=8, height=7)




# Figure A14: Conjoint results when dropping fast and slow respons --------

### subset data
d.conjoint2 <- subset(d.conjoint, d.conjoint$duration >= 0.5*median(d.conjoint$duration, na.rm=T) & d.conjoint$duration <= 2*median(d.conjoint$duration, na.rm=T))

### Regression
mm_by <- cj(d.conjoint2, 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

### prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

### define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

### get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

### plot
p <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p
ggsave("FigureA14_Pull_robustness_duration.png", p, width=8, height=7)




# Figure A15: Main conjoint results by which dimension was display --------

### Code order
d.conjoint$E2_first <- NA
d.conjoint$E2_first[d.conjoint$E2_jobs_order == 1] <- 1
d.conjoint$E2_first[d.conjoint$E2_rent_order == 1] <- 2
d.conjoint$E2_first[d.conjoint$E2_netw_order == 1] <- 3
d.conjoint$E2_first[d.conjoint$E2_immi_order == 1] <- 4
d.conjoint$E2_first[d.conjoint$E2_deport_order == 1] <- 5
d.conjoint$E2_first[d.conjoint$E2_health_order == 1] <- 6
d.conjoint$E2_first[d.conjoint$E2_type_order == 1] <- 7
d.conjoint$E2_first <- factor(d.conjoint$E2_first, levels=1:7,
                              labels = c("Jobs", "Rent", "Connections", "Share Mexican", "Deportation", "Healthcare", "Type"))


### Plot 1: MMs by order
# Regression
mm_by <- cj(d.conjoint,
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~E2_first)

# prepare output container
out <- data.frame(matrix(nrow=35*7, ncol=6))
colnames(out) <- c("first", "dimension", "attribute", "b", "lb", "ub")
out$first <- rep(1:7,each=35)
out$first <- factor(out$first, levels=1:7,
                    labels = c("Jobs", "Rent", "Connections", "Share Mexican", "Deportation", "Healthcare", "Type"))
out$dimension <- rep(rep(1:7, each=5),7)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,7)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
mm_by1 <- subset(mm_by, BY == "Jobs")
mm_by2 <- subset(mm_by, BY == "Rent")
mm_by3 <- subset(mm_by, BY == "Connections")
mm_by4 <- subset(mm_by, BY == "Share Mexican")
mm_by5 <- subset(mm_by, BY == "Deportation")
mm_by6 <- subset(mm_by, BY == "Healthcare")
mm_by7 <- subset(mm_by, BY == "Type")

# get coefficients
out[2:5,4:6] <- mm_by1[1:4,c(6,10,11)]
out[7:10,4:6] <- mm_by1[5:8,c(6,10,11)]
out[12:15,4:6] <- mm_by1[9:12,c(6,10,11)]
out[17:20,4:6] <- mm_by1[13:16,c(6,10,11)]
out[22:25,4:6] <- mm_by1[17:20,c(6,10,11)]
out[27:30,4:6] <- mm_by1[21:24,c(6,10,11)]
out[32:35,4:6] <- mm_by1[25:28,c(6,10,11)]

out[37:40,4:6] <- mm_by2[1:4,c(6,10,11)]
out[42:45,4:6] <- mm_by2[5:8,c(6,10,11)]
out[47:50,4:6] <- mm_by2[9:12,c(6,10,11)]
out[52:55,4:6] <- mm_by2[13:16,c(6,10,11)]
out[57:60,4:6] <- mm_by2[17:20,c(6,10,11)]
out[62:65,4:6] <- mm_by2[21:24,c(6,10,11)]
out[67:70,4:6] <- mm_by2[25:28,c(6,10,11)]

out[72:75,4:6] <- mm_by3[1:4,c(6,10,11)]
out[77:80,4:6] <- mm_by3[5:8,c(6,10,11)]
out[82:85,4:6] <- mm_by3[9:12,c(6,10,11)]
out[87:90,4:6] <- mm_by3[13:16,c(6,10,11)]
out[92:95,4:6] <- mm_by3[17:20,c(6,10,11)]
out[97:100,4:6] <- mm_by3[21:24,c(6,10,11)]
out[102:105,4:6] <- mm_by3[25:28,c(6,10,11)]

out[107:110,4:6] <- mm_by4[1:4,c(6,10,11)]
out[112:115,4:6] <- mm_by4[5:8,c(6,10,11)]
out[117:120,4:6] <- mm_by4[9:12,c(6,10,11)]
out[122:125,4:6] <- mm_by4[13:16,c(6,10,11)]
out[127:130,4:6] <- mm_by4[17:20,c(6,10,11)]
out[132:135,4:6] <- mm_by4[21:24,c(6,10,11)]
out[137:140,4:6] <- mm_by4[25:28,c(6,10,11)]

out[142:145,4:6] <- mm_by5[1:4,c(6,10,11)]
out[147:150,4:6] <- mm_by5[5:8,c(6,10,11)]
out[152:155,4:6] <- mm_by5[9:12,c(6,10,11)]
out[157:160,4:6] <- mm_by5[13:16,c(6,10,11)]
out[162:165,4:6] <- mm_by5[17:20,c(6,10,11)]
out[167:170,4:6] <- mm_by5[21:24,c(6,10,11)]
out[172:175,4:6] <- mm_by5[25:28,c(6,10,11)]

out[177:180,4:6] <- mm_by6[1:4,c(6,10,11)]
out[182:185,4:6] <- mm_by6[5:8,c(6,10,11)]
out[187:190,4:6] <- mm_by6[9:12,c(6,10,11)]
out[192:195,4:6] <- mm_by6[13:16,c(6,10,11)]
out[197:200,4:6] <- mm_by6[17:20,c(6,10,11)]
out[202:205,4:6] <- mm_by6[21:24,c(6,10,11)]
out[207:210,4:6] <- mm_by6[25:28,c(6,10,11)]

out[212:215,4:6] <- mm_by7[1:4,c(6,10,11)]
out[217:220,4:6] <- mm_by7[5:8,c(6,10,11)]
out[222:225,4:6] <- mm_by7[9:12,c(6,10,11)]
out[227:230,4:6] <- mm_by7[13:16,c(6,10,11)]
out[232:235,4:6] <- mm_by7[17:20,c(6,10,11)]
out[237:240,4:6] <- mm_by7[21:24,c(6,10,11)]
out[242:245,4:6] <- mm_by7[25:28,c(6,10,11)]

# plot
p <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, col=dimension), size=2,
             # shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),7),
             position=position_dodge(width=0.9)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, col=dimension), width=0,
                position=position_dodge(width=0.9)) + 
  facet_wrap(~first, ncol=4) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))
p

### save
ggsave("FigureA15_Pull_robustness_firstdimension.png", p, width=14, height=10)







# Table A12: Unweighted vs. weighted sample demographics ------------------

### declare survey design
d2 <- subset(d, !is.na(weight))
design <- svydesign(ids = ~1, data = d2, weights = ~weight)


### Get statistics
# Gender
out1 <- cbind(prop.table(table(d2$demo_gender)),c(0.5238,0.4762), prop.table(svytable(~demo_gender, design)))
out1 <- data.frame(cbind("Gender",rownames(out1),out1))
colnames(out1) <- c("variable","value","all samples", "benchmark", "weighted")

# Age group
out2 <- cbind(prop.table(table(d2$demo_agegroupALT)),c(0.1289,0.2171,0.2640,0.2049,0.1124,0.0728), prop.table(svytable(~demo_agegroupALT, design)))
out2 <- data.frame(cbind("Age group",rownames(out2),out2))
colnames(out2) <- c("variable","value","all samples", "benchmark", "weighted")

# Education
out3 <- cbind(prop.table(table(d2$demo_educationALT)),c(0.1378,0.6332,0.1133,0.0347,0.0588,0.0221), prop.table(svytable(~demo_educationALT, design)))
out3 <- data.frame(cbind("Education",rownames(out3),out3))
colnames(out3) <- c("variable","value","all samples", "benchmark", "weighted")

# Household income
out4 <- cbind(prop.table(table(d2$demo_income_quintile)),c(0.0620, 0.2220,0.2323,0.2245,0.2591), prop.table(svytable(~demo_income_quintile, design)))
out4 <- data.frame(cbind("Household income",rownames(out4),out4))
colnames(out4) <- c("variable","value","all samples", "benchmark", "weighted")

# Citizenship
out5 <- cbind(prop.table(table(d2$cit_dich)),c(0.3471,0.6529), prop.table(svytable(~cit_dich, design)))
out5 <- data.frame(cbind("Citizenship",rownames(out5),out5))
colnames(out5) <- c("variable","value","all samples", "benchmark", "weighted")


# combine all [not: Census division]
out <- rbind(out1,out2,out3,out4,out5)
rownames(out) <- NULL
out$`all samples` <- as.numeric(out$`all samples`) * 100
out$benchmark <- as.numeric(out$benchmark) * 100
out$weighted <- as.numeric(out$weighted) * 100


### Output table
print(xtable(out, caption = "Unweighted vs. weighted sample demographics", align = c("lrrccc")),
      caption.placement = "top", include.rownames=FALSE)





# Table A13: Results of the push experiment by level of education ---------

### Regressions
# baseline model
m1 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen",])

# secondary school or less
m2 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen" & 
                         d.push$demo_education %in% c("Never in school","Primary school or less", "Secondary school / GED"),])

# some college or associate's degree
m3 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen" & 
                         d.push$demo_education %in% c("Some college", "Associate's degree"),])

# College degree
m4 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen" & 
                         d.push$demo_education %in% c("Bachelor's degree", "Graduate or postgraduate degree"),])


### stargazer output
stargazer(m1,m2,m3,m4,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Pr(move away)"),
          column.labels = c("All", "HS or less", "Some college or associate's", "College degree or above"),
          covariate.labels=c("More economic opportunities",
                             "More personal connections", 
                             "Lower risk of deportation",
                             "Better access to healthcare",
                             "Fewer economic opportunities",
                             "Fewer personal connections",
                             "Higher risk of deportation",
                             "Worse access to healthcare"),
          add.lines = list(
            c("Respondent-FE?", rep("YES", 4))), 
          digits = 3)





# Figure A16: Results of the push experiment by level of education --------

### write a function (used for this and the next figure)
plotfun.appendix <- function(df){
  # create output container
  out <- data.frame(matrix(nrow=8*4, ncol=5))
  colnames(out) <- c("status","group", "condition", "b", "se")
  out$status <- rep(1:4, each=8)
  out$status <- factor(out$status, levels=1:4,
                       labels = c("all", "naturalized", "LPR", "likely unauthorized"))
  out$group <- rep(c(1,1,1,1,2,2,2,2),4)
  out$group <- factor(out$group, levels=1:2,
                      labels = c("positive\nchange", "negative\nchange"))
  out$condition <- rep(c(1,2,3,4,1,2,3,4),4)
  out$condition <- factor(out$condition, levels=1:4,
                          labels = c("Economic opportunities", "Social network", "Deportation risk", "Healthcare access"))
  
  # get values
  out[1:4,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=df))$coef[1:4,1:2]
  out[5:8,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=df))$coef[1:4,1:2]
  
  out[9:12,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=df, subset=df$cit_citizenUS_alt == "naturalized US citizen"))$coef[1:4,1:2]
  out[13:16,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=df, subset=df$cit_citizenUS_alt == "naturalized US citizen"))$coef[1:4,1:2]
  
  out[17:20,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=df, subset=df$cit_citizenUS_alt == "lawful permanent resident"))$coef[1:4,1:2]
  out[21:24,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=df, subset=df$cit_citizenUS_alt == "lawful permanent resident"))$coef[1:4,1:2]  
  
  out[25:28,4:5] <- summary(lm_robust(E1_change_bl.pos ~ I(E1_treat.pos) - 1, data=df, subset=df$cit_citizenUS_alt == "likely unauthorized immigrant"))$coef[1:4,1:2]
  out[29:32,4:5] <- summary(lm_robust(E1_change_bl.neg ~ I(E1_treat.neg) - 1, data=df, subset=df$cit_citizenUS_alt == "likely unauthorized immigrant"))$coef[1:4,1:2] 
  
  # plot
  p <- ggplot() + 
    geom_bar(data=out, aes(x=group, y=b, fill=status), 
             stat="identity", position = position_dodge(), alpha=.9) +
    geom_errorbar(data=out, aes(x=group, ymin=b-1.96*se, ymax=b+1.96*se,col=status),
                  position = position_dodge(width=0.9), width=0.25, size=1) +
    scale_fill_manual(values = c("gray50", "burlywood1","burlywood3","burlywood4"), name="") +
    scale_color_manual(values = c("gray30","gray30","gray30","gray30"), name="") +
    geom_hline(yintercept = 0, lty=2) +
    xlab("") + 
    #ylab("Change in Pr(move in next five years)\n(neg=less likely; pos=more likely)") +
    ylab("Change in Pr(move)") +
    facet_wrap(~condition) +  
    theme_bw() +
    theme(strip.text = element_text(size=14,  face="bold"),
          axis.title.x = element_blank(),
          axis.title.y = element_text(size=14, face="bold"),
          axis.text = element_text(size=12),
          legend.position = "bottom",
          legend.title = element_blank(),
          legend.text = element_text(size=12)) 
  
  # output
  return(p)  
}

### figures
# High school or less
p <- plotfun.appendix(df = subset(d, demo_education %in% c("Never in school","Primary school or less", "Secondary school / GED")))
p
ggsave("FigureA16_Push_results_educ1_HS.png", p, width=8, height=6)

# Some college
p <- plotfun.appendix(df = subset(d, demo_education %in% c("Some college", "Associate's degree")))
p
ggsave("FigureA16_Push_results_educ2_Associates.png", p, width=8, height=6)

# College degree+
p <- plotfun.appendix(df = subset(d, demo_education %in% c("Bachelor's degree", "Graduate or postgraduate degree")))
p
ggsave("FigureA16_Push_results_educ3_Collegedegree.png", p, width=8, height=6)




# Table A14: Results of the push experiment by gender ---------------------

### Regressions
# baseline model
m1 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen",])

# men
m2 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen" & 
                         d.push$demo_gender == "Male",])

# women
m3 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen" & 
                         d.push$demo_gender == "Female",])



### stargazer output
stargazer(m1,m2,m3,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Pr(move away)"),
          column.labels = c("All", "Men", "Women"),
          covariate.labels=c("More economic opportunities",
                             "More personal connections", 
                             "Lower risk of deportation",
                             "Better access to healthcare",
                             "Fewer economic opportunities",
                             "Fewer personal connections",
                             "Higher risk of deportation",
                             "Worse access to healthcare"),
          add.lines = list(
            c("Respondent-FE?", rep("YES", 3))), 
          digits = 3)



# Figure A17: Results of the push experiment by gender --------------------


### figures (uses the function created earlier)
# Men
p <- plotfun.appendix(df = subset(d, demo_gender == "Male"))
p
ggsave("FigureA17_Push_results_gender1_men.png", p, width=8, height=6)

# Women
p <- plotfun.appendix(df = subset(d, demo_gender == "Female"))
p
ggsave("FigureA17_Push_results_gender2_women.png", p, width=8, height=6)





# Figure A18: Results of the pull experiment by level of education --------

### High school or less
## Plot 1: all respondents
# regression
reg <- mm(data = d.conjoint[d.conjoint$demo_education %in% c("Never in school","Primary school or less", "Secondary school / GED"),], 
          formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- reg[1:4,c(5,9,10)]
out[7:10,3:5] <- reg[5:8,c(5,9,10)]
out[12:15,3:5] <- reg[9:12,c(5,9,10)]
out[17:20,3:5] <- reg[13:16,c(5,9,10)]
out[22:25,3:5] <- reg[17:20,c(5,9,10)]
out[27:30,3:5] <- reg[21:24,c(5,9,10)]
out[32:35,3:5] <- reg[25:28,c(5,9,10)]

# define sample
out$sample <- "all respondents"

# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2) +
  #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0) + 
  scale_color_manual(values = "chocolate4", name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("Marginal Mean") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


## Plot 2: by citizenship
# Regression
mm_by <- cj(d.conjoint[d.conjoint$demo_education %in% c("Never in school","Primary school or less", "Secondary school / GED"),], 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## Plot 3: Difference in marginal means between naturalized citizens vs. likely unauthorized immigrants
# dichotomize moderator
d.conjoint$citDich <- NA
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "naturalized US citizen"] <- 0
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "likely unauthorized immigrant"] <- 1
d.conjoint$citDich <- factor(d.conjoint$citDich, levels=0:1,
                             labels = c("naturalized", "likely unauthorized"))

# Regression
mm_diff <- cj(d.conjoint[d.conjoint$demo_education %in% c("Never in school","Primary school or less", "Secondary school / GED"),], 
              E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_differences",
              by=~citDich)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- mm_diff[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_diff[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_diff[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_diff[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_diff[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_diff[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_diff[25:28,c(6,10,11)]

# define sample
out$sample <- "difference"

# define comparison
out$comparison <- "Likely unauth. - naturalized"

# Plot
p3 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
             #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("chocolate4"), name="") +
  scale_y_continuous(limits = c(-0.25,0.25)) +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Difference in Marginal Means") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## combining all three plots
p <- ggarrange(p1,p2,p3,widths=c(0.5,0.25,0.25), nrow = 1)
p
ggsave("FigureA18_Pull_results_educ1_HS.png", p, width=13, height=7)




### Some college or Associate's degree
## Plot 1: all respondents
# regression
reg <- mm(data = d.conjoint[d.conjoint$demo_education %in% c("Some college", "Associate's degree"),], 
          formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- reg[1:4,c(5,9,10)]
out[7:10,3:5] <- reg[5:8,c(5,9,10)]
out[12:15,3:5] <- reg[9:12,c(5,9,10)]
out[17:20,3:5] <- reg[13:16,c(5,9,10)]
out[22:25,3:5] <- reg[17:20,c(5,9,10)]
out[27:30,3:5] <- reg[21:24,c(5,9,10)]
out[32:35,3:5] <- reg[25:28,c(5,9,10)]

# define sample
out$sample <- "all respondents"

# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2) +
  # shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0) + 
  scale_color_manual(values = "chocolate4", name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("Marginal Mean") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


## Plot 2: by citizenship
# Regression
mm_by <- cj(d.conjoint[d.conjoint$demo_education %in% c("Some college", "Associate's degree"),], 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## Plot 3: Difference in marginal means between naturalized citizens vs. likely unauthorized immigrants
# dichotomize moderator
d.conjoint$citDich <- NA
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "naturalized US citizen"] <- 0
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "likely unauthorized immigrant"] <- 1
d.conjoint$citDich <- factor(d.conjoint$citDich, levels=0:1,
                             labels = c("naturalized", "likely unauthorized"))

# Regression
mm_diff <- cj(d.conjoint[d.conjoint$demo_education %in% c("Some college", "Associate's degree"),], 
              E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_differences",
              by=~citDich)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- mm_diff[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_diff[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_diff[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_diff[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_diff[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_diff[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_diff[25:28,c(6,10,11)]

# define sample
out$sample <- "difference"

# define comparison
out$comparison <- "Likely unauth. - naturalized"

# Plot
p3 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
             # shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("chocolate4"), name="") +
  scale_y_continuous(limits = c(-0.25,0.25)) +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Difference in Marginal Means") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## combining all three plots
p <- ggarrange(p1,p2,p3,widths=c(0.5,0.25,0.25), nrow = 1)
p
ggsave("FigureA18_Pull_results_educ2_Associates.png", p, width=13, height=7)




### 4-year college degree or more
## Plot 1: all respondents
# regression
reg <- mm(data = d.conjoint[d.conjoint$demo_education %in% c("Bachelor's degree", "Graduate or postgraduate degree"),], 
          formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- reg[1:4,c(5,9,10)]
out[7:10,3:5] <- reg[5:8,c(5,9,10)]
out[12:15,3:5] <- reg[9:12,c(5,9,10)]
out[17:20,3:5] <- reg[13:16,c(5,9,10)]
out[22:25,3:5] <- reg[17:20,c(5,9,10)]
out[27:30,3:5] <- reg[21:24,c(5,9,10)]
out[32:35,3:5] <- reg[25:28,c(5,9,10)]

# define sample
out$sample <- "all respondents"

# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2)+
  #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0) + 
  scale_color_manual(values = "chocolate4", name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("Marginal Mean") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


## Plot 2: by citizenship
# Regression
mm_by <- cj(d.conjoint[d.conjoint$demo_education %in% c("Bachelor's degree", "Graduate or postgraduate degree"),], 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             # shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## Plot 3: Difference in marginal means between naturalized citizens vs. likely unauthorized immigrants
# dichotomize moderator
d.conjoint$citDich <- NA
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "naturalized US citizen"] <- 0
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "likely unauthorized immigrant"] <- 1
d.conjoint$citDich <- factor(d.conjoint$citDich, levels=0:1,
                             labels = c("naturalized", "likely unauthorized"))

# Regression
mm_diff <- cj(d.conjoint[d.conjoint$demo_education %in% c("Bachelor's degree", "Graduate or postgraduate degree"),], 
              E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_differences",
              by=~citDich)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- mm_diff[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_diff[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_diff[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_diff[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_diff[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_diff[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_diff[25:28,c(6,10,11)]

# define sample
out$sample <- "difference"

# define comparison
out$comparison <- "Likely unauth. - naturalized"

# Plot
p3 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
             #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("chocolate4"), name="") +
  scale_y_continuous(limits = c(-0.25,0.25)) +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Difference in Marginal Means") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## combining all three plots
p <- ggarrange(p1,p2,p3,widths=c(0.5,0.25,0.25), nrow = 1)
p
ggsave("FigureA18_Pull_results_educ3_Collegedegree.png", p, width=13, height=7)






# Figure A19: Results of the pull experiment by gender --------------------

### Men
## Plot 1: all respondents
# regression
reg <- mm(data = d.conjoint[d.conjoint$demo_gender %in% c("Male"),], 
          formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- reg[1:4,c(5,9,10)]
out[7:10,3:5] <- reg[5:8,c(5,9,10)]
out[12:15,3:5] <- reg[9:12,c(5,9,10)]
out[17:20,3:5] <- reg[13:16,c(5,9,10)]
out[22:25,3:5] <- reg[17:20,c(5,9,10)]
out[27:30,3:5] <- reg[21:24,c(5,9,10)]
out[32:35,3:5] <- reg[25:28,c(5,9,10)]

# define sample
out$sample <- "all respondents"

# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2)+
  #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0) + 
  scale_color_manual(values = "chocolate4", name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("Marginal Mean") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


## Plot 2: by citizenship
# Regression
mm_by <- cj(d.conjoint[d.conjoint$demo_gender == "Male",], 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## Plot 3: Difference in marginal means between naturalized citizens vs. likely unauthorized immigrants
# dichotomize moderator
d.conjoint$citDich <- NA
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "naturalized US citizen"] <- 0
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "likely unauthorized immigrant"] <- 1
d.conjoint$citDich <- factor(d.conjoint$citDich, levels=0:1,
                             labels = c("naturalized", "likely unauthorized"))

# Regression
mm_diff <- cj(d.conjoint[d.conjoint$demo_gender == "Male",], 
              E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_differences",
              by=~citDich)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- mm_diff[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_diff[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_diff[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_diff[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_diff[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_diff[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_diff[25:28,c(6,10,11)]

# define sample
out$sample <- "difference"

# define comparison
out$comparison <- "Likely unauth. - naturalized"

# Plot
p3 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
             #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("chocolate4"), name="") +
  scale_y_continuous(limits = c(-0.25,0.25)) +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Difference in Marginal Means") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## combining all three plots
p <- ggarrange(p1,p2,p3,widths=c(0.5,0.25,0.25), nrow = 1)
p
ggsave("FigureA19_Pull_results_gender1_men.png", p, width=13, height=7)



### Women
## Plot 1: all respondents
# regression
reg <- mm(data = d.conjoint[d.conjoint$demo_gender %in% c("Female"),], 
          formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- reg[1:4,c(5,9,10)]
out[7:10,3:5] <- reg[5:8,c(5,9,10)]
out[12:15,3:5] <- reg[9:12,c(5,9,10)]
out[17:20,3:5] <- reg[13:16,c(5,9,10)]
out[22:25,3:5] <- reg[17:20,c(5,9,10)]
out[27:30,3:5] <- reg[21:24,c(5,9,10)]
out[32:35,3:5] <- reg[25:28,c(5,9,10)]

# define sample
out$sample <- "all respondents"

# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2) +
  #shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0) + 
  scale_color_manual(values = "chocolate4", name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("Marginal Mean") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


## Plot 2: by citizenship
# Regression
mm_by <- cj(d.conjoint[d.conjoint$demo_gender == "Female",], 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm",
            by=~cit_citizenUS_alt)

# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized US citizen")
mm_by2 <- subset(mm_by, BY == "likely unauthorized immigrant")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             # shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## Plot 3: Difference in marginal means between naturalized citizens vs. likely unauthorized immigrants
# dichotomize moderator
d.conjoint$citDich <- NA
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "naturalized US citizen"] <- 0
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "likely unauthorized immigrant"] <- 1
d.conjoint$citDich <- factor(d.conjoint$citDich, levels=0:1,
                             labels = c("naturalized", "likely unauthorized"))

# Regression
mm_diff <- cj(d.conjoint[d.conjoint$demo_gender == "Female",], 
              E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_differences",
              by=~citDich)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- mm_diff[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_diff[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_diff[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_diff[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_diff[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_diff[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_diff[25:28,c(6,10,11)]

# define sample
out$sample <- "difference"

# define comparison
out$comparison <- "Likely unauth. - naturalized"

# Plot
p3 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
             # shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("chocolate4"), name="") +
  scale_y_continuous(limits = c(-0.25,0.25)) +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Difference in Marginal Means") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



## combining all three plots
p <- ggarrange(p1,p2,p3,widths=c(0.5,0.25,0.25), nrow = 1)
p
ggsave("FigureA19_Pull_results_gender2_women.png", p, width=13, height=7)



# Table A15: Results of the push experiment: weighted regression ----------

### Regressions
# baseline model
m1 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[!is.na(d.push$weight) & d.push$weight >0,], 
           weights=d.push$weight[!is.na(d.push$weight) & d.push$weight >0])

# only naturalized citizens
m2 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt == "naturalized US citizen" & !is.na(d.push$cit_citizenUS_alt) & 
                         !is.na(d.push$weight) & d.push$weight >0,], 
           weights=d.push$weight[d.push$cit_citizenUS_alt == "naturalized US citizen" & 
                                      !is.na(d.push$cit_citizenUS_alt) &
                                      !is.na(d.push$weight) & d.push$weight >0])

# only LPRs
m3 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt == "lawful permanent resident" & !is.na(d.push$cit_citizenUS_alt) &
                         !is.na(d.push$weight) & d.push$weight >0,], 
           weights=d.push$weight[d.push$cit_citizenUS_alt == "lawful permanent resident" & 
                                      !is.na(d.push$cit_citizenUS_alt) &
                                      !is.na(d.push$weight) & d.push$weight >0])

# only likely unauthorized immigrant
m4 <- felm(E1_dv ~ E1_treat | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt == "likely unauthorized immigrant" & !is.na(d.push$cit_citizenUS_alt) &
                         !is.na(d.push$weight) & d.push$weight >0,], 
           weights=d.push$weight[d.push$cit_citizenUS_alt == "likely unauthorized immigrant" & 
                                      !is.na(d.push$cit_citizenUS_alt) &
                                      !is.na(d.push$weight) & d.push$weight >0])

# interaction
m5 <- felm(E1_dv ~ E1_treat*cit_citizenUS_alt | ResponseId | 0 | ResponseId, 
           data=d.push[d.push$cit_citizenUS_alt != "natural-born US citizen" & !is.na(d.push$cit_citizenUS_alt) &
                         !is.na(d.push$weight) & d.push$weight >0,], 
           weights=d.push$weight[d.push$cit_citizenUS_alt != "natural-born US citizen" & 
                                      !is.na(d.push$cit_citizenUS_alt) &
                                      !is.na(d.push$weight) & d.push$weight >0])


### stargazer output
stargazer(m1,m2,m3,m4,m5,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = c("Pr(move away)"),
          column.labels = c("All", "Naturalized", "LPR", "Unauthorized", "All"),
          covariate.labels=c("More economic opportunities",
                             "More personal connections", 
                             "Lower risk of deportation",
                             "Better access to healthcare",
                             "Fewer economic opportunities",
                             "Fewer personal connections",
                             "Higher risk of deportation",
                             "Worse access to healthcare",
                             "LPR",
                             "Likely unauthorized",
                             "LPR X More economic opportunities",
                             "LPR X More personal connections", 
                             "LPR X Lower risk of deportation",
                             "LPR X Better access to healthcare",
                             "LPR X Fewer economic opportunities",
                             "LPR X Fewer personal connections",
                             "LPR X Higher risk of deportation",
                             "LPR X Worse access to healthcare",
                             "Likely unauthorized X More economic opportunities",
                             "Likely unauthorized X More personal connections", 
                             "Likely unauthorized X Lower risk of deportation",
                             "Likely unauthorized X Better access to healthcare",
                             "Likely unauthorized X Fewer economic opportunities",
                             "Likely unauthorized X Fewer personal connections",
                             "Likely unauthorized X Higher risk of deportation",
                             "Likely unauthorized X Worse access to healthcare"),
          add.lines = list(
            c("Respondent-FE?", rep("YES", 5))), 
          digits = 3)





# Figure A20: Determinants of destination preferences: weighted re --------

### Plot 1: all respondents
# regression
reg <- mm(data = d.conjoint[!is.na(d.conjoint$weight),], 
          formula = E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, 
          weights = ~weight,
          id = ~ResponseId)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- reg[1:4,c(5,9,10)]
out[7:10,3:5] <- reg[5:8,c(5,9,10)]
out[12:15,3:5] <- reg[9:12,c(5,9,10)]
out[17:20,3:5] <- reg[13:16,c(5,9,10)]
out[22:25,3:5] <- reg[17:20,c(5,9,10)]
out[27:30,3:5] <- reg[21:24,c(5,9,10)]
out[32:35,3:5] <- reg[25:28,c(5,9,10)]

# define sample
out$sample <- "all respondents"

# define comparison
out$comparison <- "Full sample"

# Plot
p1 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2) +
  # shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5))) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0) + 
  scale_color_manual(values = "chocolate4", name="") +
  facet_wrap(~comparison) +  
  coord_flip() + 
  ylab("Marginal Mean") +
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        #axis.text.y = element_blank(),
        axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
                                   size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))


### Plot 2: by citizenship
# dichotomize moderator
d.conjoint$citDich <- NA
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "naturalized US citizen"] <- 0
d.conjoint$citDich[d.conjoint$cit_citizenUS_alt == "likely unauthorized immigrant"] <- 1
d.conjoint$citDich <- factor(d.conjoint$citDich, levels=0:1,
                             labels = c("naturalized", "likely unauthorized"))

# Regression
mm_by <- cj(data = d.conjoint[!is.na(d.conjoint$weight),], 
            E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, 
            weights = ~weight,
            id = ~ResponseId, 
            estimate = "mm",
            by=~citDich)



# prepare output container
out <- data.frame(matrix(nrow=70, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(rep(1:7, each=5),2)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- rep(35:1,2)
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# define sample
out$sample <- rep(1:2,each=35)
out$sample <- factor(out$sample, levels=1:2,
                     labels = c("naturalized", "likely unauthorized"))

mm_by1 <- subset(mm_by, BY == "naturalized")
mm_by2 <- subset(mm_by, BY == "likely unauthorized")

# get coefficients
out[2:5,3:5] <- mm_by1[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_by1[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_by1[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_by1[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_by1[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_by1[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_by1[25:28,c(6,10,11)]

out[37:40,3:5] <- mm_by2[1:4,c(6,10,11)]
out[42:45,3:5] <- mm_by2[5:8,c(6,10,11)]
out[47:50,3:5] <- mm_by2[9:12,c(6,10,11)]
out[52:55,3:5] <- mm_by2[13:16,c(6,10,11)]
out[57:60,3:5] <- mm_by2[17:20,c(6,10,11)]
out[62:65,3:5] <- mm_by2[21:24,c(6,10,11)]
out[67:70,3:5] <- mm_by2[25:28,c(6,10,11)]

# define comparison
out$comparison <- "By immigration status"

# plot
p2 <- ggplot() + 
  geom_hline(yintercept = 0.5, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample, shape=sample), size=2,
             #shape = rep(c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),2),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("burlywood4","burlywood3"), name="") +
  scale_shape_manual(values = c(17,19), name="") +  
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Marginal Mean") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



### Plot 3: Difference in marginal means between naturalized citizens vs. likely unauthorized immigrants
# Regression
mm_diff <- cj(d.conjoint[!is.na(d.conjoint$weight),], 
              E2_choice ~ E2_jobs + E2_rent + E2_netw + E2_immi + E2_deport + E2_health + E2_type, id = ~ResponseId, estimate = "mm_differences",
              weights = ~weight,
              by=~citDich)

# prepare output container
out <- data.frame(matrix(nrow=35, ncol=5))
colnames(out) <- c("dimension", "attribute", "b", "lb", "ub")
out$dimension <- rep(1:7, each=5)
out$dimension <- factor(out$dimension, levels=1:7, 
                        labels = c("Economic opportunities", "Average rent for 1-br apt", 
                                   "Social network", "Share Mexican immigrants",
                                   "Deportation risk", "Healthcare access",
                                   "Type of community"))
out$attribute <- 35:1
out$attribute <- factor(out$attribute, levels=1:35,
                        labels = c("Neighborhood in a large city   ", "Suburban neighborhood   ", "Small town   ","Rural area   ", "Type of community",
                                   "Always   ", "Only children and pregnant women   ", "Only in emergencies   ", "Never   ", "Healthcare access for immigrants",
                                   "Protect all immigrants from deportation   ","Deport only when committed crime   ", 
                                   "No protection from deportation   ", "Trying to deport unauthorized immigrants   ", "Risk of deportation",
                                   "30%   ", "20%   ", "10%   ", "0%   ", "Mexican share of the population",
                                   "Many (10+) people   ", "Some (5-9) people   ", "A few (1-4) people   ", "Nobody   ", "Social connections",
                                   "$800   ", "$1,300   ", "$1,800   ", "$2,300   ", "Average rent 1-br apt",
                                   "Many jobs   ", "Some jobs   ", "Very few jobs   ", "No jobs   ", "Job opportunities"))

# get coefficients
out[2:5,3:5] <- mm_diff[1:4,c(6,10,11)]
out[7:10,3:5] <- mm_diff[5:8,c(6,10,11)]
out[12:15,3:5] <- mm_diff[9:12,c(6,10,11)]
out[17:20,3:5] <- mm_diff[13:16,c(6,10,11)]
out[22:25,3:5] <- mm_diff[17:20,c(6,10,11)]
out[27:30,3:5] <- mm_diff[21:24,c(6,10,11)]
out[32:35,3:5] <- mm_diff[25:28,c(6,10,11)]

# define sample
out$sample <- "difference"

# define comparison
out$comparison <- "Likely unauth. - naturalized"

# Plot
p3 <- ggplot() + 
  geom_hline(yintercept = 0, lty=2) +
  geom_point(data=out, aes(x=attribute, y=b, color=sample), size=2,
             # shape = c(rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5), rep(17,5), rep(19,5)),
             position=position_dodge(width=0.8)) +
  geom_errorbar(data=out, aes(x=attribute, ymin=lb, ymax=ub, color=sample), width=0,
                position=position_dodge(width=0.8)) + 
  scale_color_manual(values = c("chocolate4"), name="") +
  scale_y_continuous(limits = c(-0.2,0.2)) +
  facet_wrap(~comparison) +
  coord_flip() + 
  ylab("Difference in Marginal Means") + 
  theme_bw() +
  theme(strip.text = element_text(size=14, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=11),
        axis.text.y = element_blank(),
        #axis.text.y = element_text(face = rep(c(rep("plain",4),"bold"),6),
        #                           size = rep(c(rep(11,4),11),6)),
        axis.ticks.length.x  = unit(0.2, "cm"),
        axis.ticks.length.y  = rep(unit(0.2, "cm"),4))



### combining all three plots
p <- ggarrange(p1,p2,p3,widths=c(0.5,0.25,0.25), nrow = 1)
p
ggsave("FigureA20_Pull_results_weighted.png", p, width=13, height=7)











